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ABSTRACT 

Within  the  past  decade  easy  access  to  high  speed  digital  computers  has 
renewed  interest  in  deriving  molecular  weight  distributions  from  sedimentation- 
diffusion  equilibrium  data.  One  of  the  computational  schemes  which  appears 
most  promising  is  the  Simplex  Method  of  linear  programming.  The  purpose  of 
this  work  was  to  investigate  the  advantages  and  limitations  of  this  approach. 

It  was  found  that,  even  though  inferring  a  molecular  weight  distribution 
from  sedimentation-diffusion  equilibrium  data  is  mathematically  an  ill-posed 
problem,  the  method  of  linear  programming  yields  qualitatively  a  good  molec¬ 
ular  weight  distribution.  Also,  the  method  proved  satisfactory  for  the  case 
when  sedimentation  equilibrium  data  was  acquired  from  only  a  single  angular 
velocity. 
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SECTION  I 
INTRODUCTION 

The  relationship  describing  sedimentation-diffusion  equilibrium  of  an  ideal 
polydisperse  solution  in  an  ultracentrifuge  can  be  given  by  a  Fredholm  integral 
equation  of  the  first  kind  (Reference  1).  Since  no  rigorous  solution  of  this  in¬ 
tegral  is  known,  there  have  been  many  attempts  to  solve  it  by  approximation 
(References  2,  3).  These  efforts  mainly  involved  use  of  Fourier  transforms  or 
Laplace  transforms  by  assuming  an  approximate  functional  expression  for  the 
experimental  concentration  gradient  along  the  ultracentrifugal  cell,  or  by  ex¬ 
panding  the  molecular  weight  distribution  (MWD)  into  a  polynomial  of  assumed 
functions. 

The  main  weakness  has  been  that  some  parts  of  the  calculated  distribution 
would  be  negative.  Physically,  of  course,  we  know  that  the  MWD  for  any 
molecular  weight  must  always  be  positive  or  zero.  Recently  Lee  (Reference  4) 
carried  out  an  investigation  of  the  Fredholm  integral  equation  and  found  that 
mathematically  it  is  an  "ill-posed"  problem.  In  trying  to  infer  a  MWD  from 
experimental  measurements  of  concentration  gradients  small  errors  can  lead 
to  an  unacceptable  MWD.  Therefore,  we  compromised  in  trying  to  determine 
only  an  "overall"  shape  of  the  MWD  without  being  specific  to  individual  points, 
i.  e. ,  we  allowed  certain  fluctuations  of  the  curve  to  be  present  and  ignored 
fine  structure. 

To  generalize  a  theoretical  analysis,  let  us  accept  that  the  MWD  can  be 
slightly  negative  for  some  molecular  weight  values.  Since  we  chose  to  ignore 
the  point-by -point  functional  form  of  the  MWD,  the  next  logical  step  would  be 
to  subdivide  the  MWD  into  narrow  (not  infinitesimal)  but  finite  molecular  weight 
strips.  This  would  result  in  approximations  of  MWD  by  rectangles  of  finite 
width  and  would  lead  naturally  to  the  use  of  matrices.  This  has  been  done 
(Reference  5)  but  unfortunately  the  matrices  are  "ill-conditioned"  or  nearly 
singular. 

Scholte  (Reference  6  and  7)  in  1968,  still  using  matrices,  applied  the 
scheme  of  linear  programming  to  infer  a  MWD  from  experimental  measure¬ 
ments  of  concentration  gradients  at  various  angular  velocities.  The  main 
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advantage  to  this  approach  is  that  values  of  the  MWD  are  forced  to  be  greater 
than  or  equal  to  zero  and  "slack  variables"  are  introduced  to  account  for  ex¬ 
perimental  error.  Scholte  evaluated  the  MWD  at  ten  molecular  weights,  then 
shifted  to  ten  other  moLeeular  weights  in  a  prescribed  manner,  continuing 
until,  finally,  there  were  four  such  sets.  Since  each  set  represented  an  indi¬ 
vidual  solution,  one  quarter  of  the,  sum  of  the  four  sets  also  represented  a 
solution.  By  doing  this,  Scholte  obtained  good  agreement  between  his  assumed 
and  calculated  molecular  weight  distributions. 

There  are,  however,  three  reasons  why  Scholte's  scheme  cannot  be  blindly 

applied  to  other  systems.  These  are:  (1)  Scholte  dealt  with  a  molecular  weight 
4  6 

range  of  5  x  10  to  10  ;  by  comparison,  in  many  cases  of  synthetic  polymers 

5 

the  range  is  much  narrower,  e.  g. ,  0  to  10  .  (2)  Scholte  used  five  or  more 
angular  velocities,  each  requiring  several  days  for  equilibrium.  There  are, 
however,  cases  when  equilibrium  at  each  velocity  requires  a  much  longer 
time  (Reference  8).  Therefore,  it  is  important  to  have  a  scheme  which  would 
produce  a  MWD  from  data  taken  at  one  velocity.  (3)  Since  our  interest  was  in 
a  different  molecular  weight  range  and  we  were  using  experimental  data  from 
only  one  angular  velocity,  the  effects  of  experimental  error  on  the  calculated 
MWD  had  to  be  investigated. 

A  computer  program  using  Scholte’s  ideas  was  independently  coded  and  a 
different  linear  programming  (LP)  solving  routine  was  employed.  The  new 
program  reliability  was  verified  by  reproducing  Scholte's  published  results. 
Then  application  of  the  new  program  to  new  specific  needs  stated  above  were 
investigated.  A  brief  description  of  linear  programming  theory  follows. 
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SECTION  II 
THEORY 


In  the  brief  discussion  which  follows,  all  theorems  and  definitions  are  given 
without  proof  or  examples.  All  material  on  linear  programming  was  taken  from 
other  works  (References  9,  10,  and  11). 

Definition  1.  A  simplex  is  an  n-dimensional  convex  polyhedron  having 
exactly  n+1  vertices.  The  boundary  of  the  simplex  contains  simplical  faces 
of  dimension  i  where  i  <  n.  The  number  of  such  faces  of  dimension  i  is  ( .+j) 
where  (^)  =  n!/m!(n-m)!.  A  simplex  in  zero  dimension  is  a  point,  in  1- 
dimension  a  line,  in  2-dimension  a  triangle,  in  3-dimension  a  tetrahedron, 

etc.  The  equation  of  a  simplex  with  unit  intercept  is  X.  >  0  and  Z  X.  =  1. 

1  i  1 

Definition  2.  A  subset  C  of  (n-dimensional  Euclidean  Space)  is  a  convex 
set  if  and  only  if  for  all  pairs  of  points  and  V2  in  C  any  convex  combination 

v  =  /3,  v,  ♦  0 2v,  (l) 

is  also  in  C,  where  /3  .  are  scalars,  /3  >0,  and  Z  =  1. 

X  1  j  1 

Definition  3.  A  point  V  in  a  convex  set  C  is  called  an  extreme  point  if  V_ 
cannot  be  expressed  as  a  convex  combination  of  any  other  two  distinct  points  in 
C.  That  is,  if  we  denote  the  convex  set  of  solutions  to  the  linear  programming 
problem  by  K  and  if  K  is  a  convex  polygon,  then  K  is  the  convex  hull  of  the 
extreme  points  of  K.  Therefore,  every  feasible  solution  in  K  can  be  represented 
as  a  convex  combination  of  the  extreme  feasible  solutions  in  K. 

Theorem  1.  The  set  of  all  feasible  solutions  to  the  linear  programming 
problem  is  a  convex  set. 

In  general,  the  linear  programming  problem  can  be  described  as  follows: 
Given  is  a  convex  set  defined  by  a  set  of  linear  constraints  in  E  .  From  all  the 
points  belonging  to  the  convex  set,  we  wish  to  determine  a  subset  of  points 
(which  will  contain  either  one  or  many  points)  for  which  a  linear  objective 
function  is  optimized. 
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Usually  we  are  confronted  with  a  set  of  simultaneous  equations 


a  IIXJ  +  °I2X2+  °I3X3  + 


+  a  i  m  x  m  =  b  i 


ani  x ,  +  a__x+a__x_  + 
ni  i  m  2  n3  3 


+  anmxm  “  bn 


(2) 


where  n  >  m.  For  simplicity  let  m  equal  n;  let  A  be  the  matrix,  j  ,  (i  = 
1, 2, 3,  — mand  j  =  1,  2, 3 — m);_X  the  vector  [x^  j  ,  j  =  1,  2,  3 — m  and  b  the 
vector  |  b.  j  ,  i  =  1 , 2,  3 — ,  m.  Then  Equation  2  can  be  written  in  the  form 


Ax  =  b 


(3) 


Since  A  is  a  square  matrix  and  assumed  nonsingular,  the  solution  vector  is 
expressed  as 

_x  =  A-1b  (4) 

A  simple  computational  scheme  is  the  complete  elimination  method  of  Jordan 
and  Gauss  which  has  a  finite  number  of  steps  or  iterations.  In  just  m  iterations 
the  procedure  multiplies  the  system  (Equation  2)  by  A  1  to  obtain  Equation  4. 
This  is  the  standard  matrix  problem  which  is  assumed  familiar  to  the  reader. 

Now,  in  linear  programming  the  problem  is  reversed.  Instead  of  having 
an  "over-determined"  system  as  indicated  by  Equation  2,  we  have  an  "under¬ 
determined"  system  (i.  e. ,  n  <  m)  subject  to  other  constraints.  That  is,  we 
wish  to  find  a  vector  j  x^  j  ,  i  =  1, 2,  3,  — m  which  minimizes  the  linear  form 
(i.  e. ,  the  objective  function) 


C1X1  +  C2X2  +  C3X3+  " 


-  +  C  x 
m  m 


subject  to  linear  constraints 

x.  >  o,  j  =  1,2,3,—  ,m 


(5) 

(6) 


and  the  set  of  equations  given  by  Equation  2  but  with  n  <  m. 


For  large  n  and  m  it  would  be  an  impossible  task  to  evaluate  all  possible 
solutions  and  select  one  that  minimizes  the  objective  function.  A  computational 
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scheme  is  desired  which  converges  to  a  minimum  solution.  The  Simplex 
Method,  devised  by  Dantzig  (Reference  11)  is  such  a  scheme.  In  Reference  11  the 

equation  Z  Xi  is  used  as  a  constraint.  The  procedure  finds  an  extreme  point 
and  determines  whether  it  is  the  minimum.  If  it  is  not,  the  procedure  finds 
a  neighboring  extreme  point  whose  corresponding  value  of  the  objective  function 
is  less  than  or  equal  to  the  preceding  value.  In  a  finite  number  of  such  steps 
(usually  between  n  and  2n),  a  minimizing  feasible  solution  is  found.  The  Sim¬ 
plex  Method  makes  it  possible  to  discover  whether  the  problem  has  any  finite 
minimizing  solutions  or  no  feasible  solutions  at  all. 

Consideration  is  now  given  to  how  this  can  be  related  to  the  problem  at 
hand,  namely,  molecular  weight  determination  via  sedimentation-diffusion 
equilibrium.  The  equation  describing  sedimentation-diffusion  equilibrium  for 
a  heterogeneous  system  is  given  (from  Reference  1)  by 


where. 


(7) 


(8) 


In  the  above  equations  C°  is  the  concentration  of  the  original  solution,  C  is  the 

equilibrium  concentration  at  radial  distance  r,  M  is  the  molecular  weight, 

F(M)  is.  the  frequency  function  of  molecular  weight,  £  =  (r^  -r  )/(r^  -rffl  ) 

with  r  the  radial  distance  from  the  center  of  rotation  to  the  meniscus,  and  r, 
m  ’  b 

the  radial  distance  from  the  center  of  rotation  to  the  bottom  of  the  cell.  Also, 

2  2  2 

X  =  (1-v/d  )  a)  (r^  -r  )/2RT,  where  v  is  the  partial  specific  volume  of  the 
dissolved  substance,  p  is  the  density  of  the  solution,  u)  is  the  angular  velocity 
in  radian  per  second,  R  is  the  gas  constant,  and  T  the  absolute  temperature. 


Rewriting  Equations  7  and  8  for  the  discrete  case  (Dirac  8  -functions)  one 
obtains 


u(  M„>  =  l 


Xf  MS,  e 


-  *■! 

i  m^n 


I  _  e~  Mm 


•m 


(9) 


5 


AFML-TR-67-121 
PART  III 

and, 


(10) 


where 


U(Alf£j 


d_ 


and 


f 


m 


is  the  weight  fraction  of  molecules  of  a  given  molecular  weight  in  the 
original  sample.  Recall  that  the  U  (  Xj,  £n)  and  £n  are  the  experimentally 
measured  quantities  with  X .  being  the  product  of  a  constant  (determined  from 
auxiliary  measurements)  and  the  square  of  the  angular  speed  of  the  rotor. 


For  convenience  of  notation  let 


K 


.2  2  -X.  M  £ 

X,  M_  e  i 


i  '"m 


£n 


i-  rxi 


and 


ul  *  u  <  Xi,«„  1 


(I  I) 


(12) 


where  for  each  i,  n  =  1, 2,  — N;  i  =  1,2,  — ,1;  m  =  1, 2,  — M;  and  ■£=  1, 2,  — , 
L  with  L  =  IN  and  L  >  M. 


Thus  Equation  9  becomes 


U. 


=  I  K 


m 


4m  *m 


(13) 


Since  the  quantities  of  Ug  are  ejqjerimentally  measured,  they  will  in  all  prob¬ 
ability  be  greater  than  or  less  than  their  true  precise  value  (i.e.,  there  exists 
experimental  error).  Although  this  physical  fact  is  accepted,  experimentally 
Equation  13  does  not  hold  true.  This  is  especially  apparent  when  we  investigate  the 
matrix  j  Kgmj  and  find  it  ill-conditioned.  In  essence,  the  matrix 
acts  as  an  "amplifier"  for  any  error  which  might  exist  in  the  set 


If  we  grant  that  an  error  in  Ug  exists,  Equation  13  becomes 


=  I 

m 


K  -  f  + 
«m  m 


(14) 
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where  is  the  experimental  error  in  .  Since  the  application  of  linear  pro¬ 
gramming  necessitates  that  all  (see  Equation  1)  are  positive  or  zero,  we  must 
account  for  error’s  being  positive  or  negative.  It  is  the  inclusion  of  error  that 
now  enables  us  to  go  from  an  ’’over  determined"  system  to  an  ’’under  deter¬ 
mined"  system.  The  linear  programming  procedure  is  now  applicable.  In 
particular,  you  will  recall  that  a  modified  Simplex  Method  can  be  used. 

Recapitulating,  we  now  obtain  the  formulation  of  the  linear  programming 
scheme  as  used  to  determine  the  MWD  from  sedimentation-diffusion  equilib¬ 
rium.  We  wish  to  find  the  set  j  fm  j  ,  m  =  1, 2,  — Q,  which  minimizes  the 
linear  form  (i.  e. ,  the  objective  function) 

L 

I  (8.  +  fig)  (15) 

g-\  JL 


subject  to  the  linear  constraints 


and 


fm  *  0, 


m 


s  °}J 
>  0  J 


1.2,3, 

1.2,3, 


•  Q 

•  L 


(16) 


V.  +  K.2f2+  *  •  +  K,QfQ+  S.-£.  +  0  +  0+-  •  +  0  +  0  =  U, 

K2,f,  +  K22f2+  •  •  •  +  K2QfQ+  0  +  0  +  S2-  &z+  ■  •  •  +0  +  0  =  U2 

!  (17) 

i 

i 

V/W  +klo’o+  o  +  o+o+o+--  +Sc0L=uL 

where  L  >  Q**.  Let  the  set  j  x.  j  ,  i  =  1, 2, 3, . .  .,Q  +  2L  be  composed  of  f 

values  for  i  =  1,2, 3,  — Q,  and  alternately  8^  and  ftgior  i=Q  +  l,  Q  +  2, - , 

Q  +  2L;  where  all  >  O.  Also  let  the  L  x(Q  +  2L)  matrix  P  be  represented  by 


**This  is  not  an  absolutely  necessary  condition.  When  one  velocity  was  used 
the  situation  arose  where  L  <  Q. 
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K 

K  *  * 

*  K  A 

il 

12 

|Q 

K 

K  •  * 

•  K  « 

21 

22 

\ 

20 

1 

\ 

1 

1 

\ 

1 

X 

X 

1 

X 

X 

Ku 

K.  * 

Lz 

•  •  klq 

I 

0 

i 

i 

i 

i 

i 

0 


Then,  In  matrix  notation  the  problem  is  formulated  by 

PX  =  U 

The  next  section  will  describe  the  application  of  this  method. 


(18) 


(19) 


8 


AFML-TR-67-121 
PART  HI 

section  m 

VARIABLE  FACTORS  IN  COMPUTATION  AND 
THEIR  INFLUENCE  ON  RESULTING  MWD 


A.  FORMULATION  OF  COMPUTER  PROBLEM 


The  objective  of  this  section  is  to  present  the  results  which  three  variable 
factors  investigated  have  on  a  MWD.  Since  the  actual  programming  involved  a 
slight  modification  of  Equation  9,  a  listing  and  discussion  of  the  variable 
factors  investigated  will  be  preceded  by  a  discussion  of  the  actual  equation 
programmed. 


The  programmed  equation  is  given  by 


u(xk,e.i>  ■  i 


X2  M  2  e  “Xr  M"  f(m„) 


n-  AMn 


_  XkMn 


(20) 


I  -  e 


Here  all  quantities  have  the  same  meaning  as  in  Equation  7  and  9.  However, 
we  must  remember  that,  since  X  is  proportional  to  the  square  of  the  angular 
velocity,  the  index  k  indicates  the  various  velocities  at  which  equilibrium  was 
achieved.  For  each  velocity  there  exists  a  set  of  £ -values,  i.  e. ,  for  each  X  ^ 
there  is  a  corresponding  set  {Cj}  •  If  there  are  data  from  five  velocities  and 
for  each  velocity  there  corresponds  five  £  -  values ,  this  would  imply  twenty- 
five  U  -values.  It  is  imperative  that  the  molecular  weight  range  being  investi¬ 
gated  incorporate  all  molecular  weights  present  in  the  solution  sample  being 
centrifugated.  Since  Equation  9  and  20  deal  with  discrete  molecular  weights, 
some  procedure  must  be  employed  to  span  the  entire  molecular  weight  range 
(MWR).  FoUowing  the  idea  of  Scholte  (References  6  and  7),  a  multiplicative 
factor  (g -factor)  has  been  introduced.  Therefore,  starting  with  the  first 
molecular  weight  M^,  the  other  molecular  weights  in  a  given  sampled  set 
could  be  generated.  Knowing,  a  priori,  the  MWR  we  can  now  calculate  the  g- 
factor  and  the  number  (N„)  of  molecular  weights  needed  to  span  the  MWR. 

(The  g-factor  is  later  used  as  a  variable  parameter  related  to  the  error  in 
the  experimental  concentration  gradients.  However,  the  value  of  g  must  be  at 
least  large  enough  to  ensure  the  actual  MWR  present  in  a  given  experiment 
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is  spanned.  A  convenient  technique  for  finding  the  MWR  is  given  in  Reference 
8).  The  molecular  weights  sampled  will  be 

M  =  M  gn"'  (2D 

n  i  "  v  ' 


where  n  =  1, 2,  3,  — ,Nq.  This  enables  one  to  divide  the  MWR  into  non-over¬ 
lapping  subranges.  Each  subrange  span  is  denoted  by 


AM  =  M  - M 
n  n  n  -i 


(22) 


where  Mq  is  assumed  zero.  By  employing  the  averaging  technique  of  Scholte 
(Reference  7),  after  solving  Equation  20  for  one  set  of  molecular  weights, 
a  new  set  of  molecular  weights  is  selected  in  a  prescribed  manner.  The  new 
set  is  shifted  relative  to  the  previous  set  by  a  multiplicative  factor  g  where 
N  is  the  number  of  desired  molecular  weight  sets.  That  is,  if  the  number 
of  the  set  is  labeled  by  the  index  n  and  the  molecular  weights  within  a  set  by  j, 
then 

Mj,n  *  <23> 

and 

AM.  =  (  g'_  g',/2)  M.  (24) 

jn  jn 


where  j  =  1,2,3, - ,NQ;  n  =  1,2,3, - ,N;  and  Mq  n  =  O. 

The  concentration  of  molecular  weights  in  a  given  subrange  is  simply  the 

weight  fraction  (f  )  multiplied  by  the  initial  solution  concentration  (C°  =  C°f  ). 

m  m  nr 

When  Scholte  (Reference  7)  presents  his  final  results  they  are  in  the  form  MF(M) 
versus  M.  In  this  work  a  modified  system  F(M)  versus  M  has  been  calculated 
(Equation  20),  since  one  usually  has  less  qualitative  feeling  for  MF(M)  than  for 
F(M). 

Now,  it  is  possible  to  list  the  variable  parameters  investigated  in  this 
work.  They  are  as  follows: 

1.  The  effect  that  varying  the  g-factor  (i.  e. ,  the  span  of  the  molecular 
weight  subrange)  has  on  the  resulting  MWD;  and  also  the  effect  when  the  number 
of  sets  of  molecular  weights  sampled  was  varied. 
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2.  How  the  resulting  MWD  is  affected  by  varying  the  number  of  elements 
in  the  sets  {£.  }  and  |\jJ  ’ 

3.  Whether  the  introduction  of  error  into  idealized  U-values  affects  the 
resulting  MWD.  This  includes  normal  random  error  and  weighted  random 
error. 

The  results  of  each  will  be  discussed  successively  in  the  following 
section. 

B.  DISTRIBUTIONS  STUDIED 

The  g-factor  is  related  to  the  experimental  error  in  the  u-values  (Ref¬ 
erence  7).  Figures  1  through  4  show  the  effect  of  varying  the  g-factor  and  the 
number  of  sets  of  molecular  weights.  In  each  case  the  solid-line  curve  is  the 
assumed  distribution  from  which  u-values  were  calculated.  For  all  curves, 
seven  velocities  were  used  and  with  each  velocity  five  £  -values.  The  squared 
angular  velocities  were  4. 1693  x  10^,  5.  8370  x  105,  8. 1718  x  105,  11. 4406  x 
105,  16.  0168  x  105,  22.4235  x  105,  and  31.  3929  x  105  rad2/sec2  with  £  =  O, 
1/4,  12,  3/4,  and  1  for  each  case.  Therefore,  thirty-five  U-values  were  used 
for  these  calculations. 

As  previously  mentioned,  we  must  always  be  sure  that  the  MWR  is  wide 
enough  to  incorporate  all  molecular  weights  present  in  the  sample.  To  study 
the  effect  of  range  size  on  the  MWD,  all  parameters  in  Figure  4  were  held 
constant  except  the  value  for  MWR,  which  was  brought  closer  to  the  MWR  of 
the  assumed  distribution.  As  shown  in  Figure  5,  structure  begins  to  appear 
when  the  highest  molecular  weight  sampled  was  not  far  beyond  the  actual 
highest  molecular  weight  present.  Figure  5  represents  a  calculation  involving 
twenty  molecular  weight  sets.  When  the  number  of  molecular  weight  sets  was 
decreased  from  twenty  to  ten  (Figure  6)  then  five  (Figure  7),  there  was  no 
appreciable  change  except  that,  naturally,  the  calculated  points  were  spaced 
farther  apart. 

What  would  be  the  effect  if  the  number  of  £  -values  associated  with  each 
velocity  was  increased?  As  previously  mentioned,  five  £  -values  have  been 
used  per  velocity  for  Figures  1  through  7.  Figure  8  shows  the  results  with  all 
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parameters  of  Figure  6  held  constant  except  that,  now,  nine  £  -values  were 
used.  The  nine  £  -values  were  so  chosen  that  £  =  O  to  1  with  A  £  =  1/8.  From 
the  study  of  this  assumed  distribution  it  appears  that  using  five  £ -values, 
sampling  twenty  sets  of  molecular  weights,  and  using  a  g-f actor  ~  2.  0  seemed 
to  have  produced  the  optimum  desired  results,  i.  e. ,  the  calculated  MWD 
agreeing  best  with  the  assumed  MWD.  If  the  g-factor  became  too  small  the 
result  was  noise,  that  is,  the  accuracy  of  the  U-values  did  not  warrant  such 
precision,  or  the  matrix  in  the  LP  solver  routine  became  singular. 

At  this  stage  a  normal  Gaussian  distribution  with  a  MWR  of  0  to  120, 000 
was  investigated.  Once  again  the  g-factor  was  varied.  The  lowest  g-value  used 
was  1. 15  and  the  highest  4. 0.  The  former  resulted  in  an  erratic  MWD  and  the 
latter  resulted  in  a  curve  which  went  exponentially  to  zero  at  high  molecular 
weights.  The  best  g-value  for  this  specific  case  was  g  =  1. 8.  In  general,  a 
satisfactory  technique  was  to  start  with  a  low  value  of  g.  Then  as  the  value  of 
g  was  increased  the  erratic  behavior  of  the  MWD  disappeared.  At  the  g-value 
where  the  erratic  results  seemed  to  disappear,  that  value  was  established  as 
the  appropriate  one.  Then  the  maximum  reliable  "fine  structure"  for  a  given 
set  of  experimental  U  -values  was  attained. 

Using  this  assumed  normal  distribution  (its  functional  form)  the  MWR  was 

shifted  to  investigate  the  reliability  of  the  method  for  various  molecular  weight 

5  6 

ranges.  One  range  tried  was  from  0  to  12, 000  and  another  from  10  to  10.  In 
each  case  the  results  were  satisfactory,  considering  that  in  all  cases  g  was 
kept  constant  (g  =  1.  8). 

As  previously  mentioned,  the  u-values  used  resulted  from  seven  assumed 
velocities  ranging  from  about  6, 000  to  50, 000  RPM.  It  would  be  advantageous 
if  the  experimental  U's  were  obtained  from  an  equilibrium  sedimentation- 
diffusion  experiment  at  only  one  angular  velocity.  To  check  this,  the  normal 
Gaussian  distribution,  MWR  from  0  to  120,000,  was  approximated  (holding  all 
other  variables  constant)  by  deleting  all  U-values  associated  with  various 
angular  velocities.  All  combinations  of  the  velocities  were  tried.  By  using 
only  the  lowest  angular  velocity  (6, 166  RPM),  the  computer  program  produced 
a  MWD  which  "fit"  the  assumed  MWD  as  well  as  the  case  where  all  seven 
angular  velocities  were  used.  In  fact,  all  single  velocity  cases  resulted  in 
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a  reasonable  MWD.  Therefore,  we  would  conclude  that  an  acceptable  MWD 
could  be  obtained  from  an  equilibrium  sedimentation-diffusion  experiment 
at  one  angular  velocity,  at  least  with  a  MWR  of  0  to  120, 000. 

The  third  area  of  investigation  involved  use  of  the  normal  Gaussian  MWD 
(0  <  M  <  120, 000)  U- values  from  one  angular  velocity  (6, 166  RPM)  and  at 
£  =  0,  1/4,  1/2,  3/4,  and  1  to  find  the  effect  that  error  in  the  U-values  would 
have  on  the  calculated  MWD.  Error  (1,  2,  5,  10,  and  20%,  respectively)  was 
introduced  by  aid  of  a  random  error  generator.  For  each  magnitude  of  error, 
five  calculations  were  employed  to  vitiate  any  wrong  conclusions  that  one  error 
distribution  might  have  on  the  final  MWD.  For  each  case  (1  to  20%)  when  the 
error  was  normally  distributed  (dotted  line  Figure  9)  the  calculated  MWD 
agreed  with  the  assumed  MWD.  Naturally  the  1%  error  case  gave  the  best 
"fit"  to  the  assumed  MWD,  but  even  for  the  20%  error  case  the  calculated 
MWD  was  not  unacceptable.  When  the  error  introduced  in  the  U  -values  was 
such  as  to  be  weighted  (dashed  curves  Figure  9),  the  calculated  MWD  was 
entirely  different  from  the  expected  normal  MWD.  This  phenomenon  sub¬ 
stantiates  the  findings  of  Lee  (Reference  5)  and  Tikhonov  and  Glasko  (Ref¬ 
erence  12). 
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SECTION  IV 

DESCRIPTION  OF  COMPUTER  PROGRAM 
FLOW  DIAGRAM 
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DESCRIPTION  OF  VARIABLES 


a.^  -  entry  in  matrix  for  LP  problem 

j 


m 


b. 

1 

C° 

c° 

f 

n 


F. 

g 


n 

H 

m 

M 

n 

M. 

N 

Nh 

N 

Ncu 

N 

V 

xj <L 
Xk 

Ah 

Am. 

J.n 


cu 


-  largest  x  value 

-  right-hand  side  of  LP  problem 

-  initial  concentration  of  solution 

h 

-  initial  concentration  of  solution  whose  molecular  weight  is -77- 

-  weight  fraction  of  molecular  weight 

-  results  of  LP  solution,  frequency  list  for  mol.  wt.  of  ^ 

-  M  multiplier 

-  input  test  array  which  is  a  function  of  molecular  weight 
(References  5  and  8) 

=  2RT 
1  -vp 

-  smallest  x  value 

-  molecular  weight  array  associated  with  input  test  array 

-  molecular  weight  matrix  of  values  used  in  LP  solution 

-  number  of  LP  sets  to  try 

-  number  of  input  h  values 

-  number  of  M  values  to  use  for  use  LP  solution 

2 

-  number  of  input  uj  values 

-  number  of  input  x  values 

1  dc 
“  C°  d  £ 

-  array  of  distances  squared  from  center  of  rotation 

2 

-  function  of  ai  k 

-  constant  difference  between  values  of  h^  array 

-  difference  between  successive  i  values  of  M.  matrix 

J 

-  function  of  x  ^ 

-  square  of  the  angular  velocity 
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IN  1  ROUTINE 

Input:  Nx,  h^,  Ah,  Nn>  H,  C°,  gj^2  array,  where  k  =  1, 

I 

x.  array ,  where  Si-  1, 
C°n  array,  where  n  =  1, 


X*  hn+l  =  hn  +  Ah  for  n  =  !» — ,  (Nh-1) 
2.  m  =  x^ 


X  k  =  (b-m)  •  H  •  wk2  for  k  =  1, — ,NW 


3. 

4. 


b-x^ 
^  Si  ~  b-m 


iovl=  1, - ,Nx 


h 

if”  for  n  =  1, — ,Nh 


C° 

f  =  —a. 

n  C° 


U 


M 


Nh 

=  I 

n  =  I 


(  Xk  Mn  f  e  XkMn 

'  X,,  M 


-  e 


'k  n 


6. 


for  k  =  1, — ,NW  and  1=1, — ,Nx 
2 

Write  out  ai  array,  £  array,  and  U  matrix 


7.  Call  LP  SOLVER 
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IN  2  ROUTINE 

2 

Input:  N.  ,  N  ,  H,  C°,  oj*  array  where  k  =  1, — ,N 

^  X  K  U) 

x,  array  where  4,  =  1,  — ,N 
it  x 

dC 

.  matrix  where  k  -  1,  — ,N,j  and 
fNx 


2. 


b  = 

X  k 


n 


x 

=  (b-m) 


for  k  = 


.  b"xA 

3-  ^  =b-m  £orA=1--“Nx 

4-  uk,*=C°  Wk,t  for  k  =  1,  ,NW 

and  1, - ,N 

X 

2 

5.  Write  out  u>  array,  £  array,  and  U  matrix 

6.  Call  LP  SOLVER 
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LP  SOLVER 

1.  Input  from  IN  1:  N^,  N  ,  X  array,  £  ai'ray,  and  U  matrix 

2.  Input  from  cards:  M.]^  1,  g,  N^,  N,  HEADING 

3.  Write  out  input  from  cards 


4.  N  N  r  N 
row  OJ  x 


N  ,  =  N_.  +  2  N 
col  M  row 

n  =  1 

5.  Calculate  the  following  matrix  entries 

a.  =  1  for  i  =  1,  —  ,Nrow  and  j  =  NM  +  1,  Ncol_1>  2 

j 

a.  =  1  for  i  =  1,  — ,Nrow  and  j  =  NM  +  1,  g 

_  <XkMj,n>2 


a.  = 

i. 

J 


AM 


1  -  e  k  3,n 


for  i  =  1, — »Nrow  and  j  =  1, - ,N. 


M 


where  k  =  1  |  +  1 

t=  i  (k-1), •  N 


-1 


and 


M.  =  M-.  .  g^ 

j,n  l,n 


AM.in  =  fe1/2-g-1/2)Mjin 


6.  Calculate  the  following  right-hand  sides  b.  =11^  where  i,k,  and  SL 
are  defined  as  in  step  5. 

7.  Calculate  the  following  objective  coeffcients 

C.  -  O  for  j  =  1, - ,NM 

O.  =  1  for  J=NM ■+!,-- ,Nool 
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LP  SOLVER 


8.  Write  out  the  determinants  of: 

[a. J  for  i  =  1, — ,Nm  and  j  =  1, — tNj 

[a.Jfori  =  NM  +  l, — ,2-  NM 


a.  for  aH  sets  of  N, ,  rows  less  than  N 
1  i.J  M  row 

J 

9.  Write  out  matrix  fa.  1  for  i  =  1,  — >Nrow  aad  j  =  1,  — ,N^  in  exponent 

form.  j J 

10.  CaULP  solver  and  store  solutions  in  F.  array. 

J,n 

11.  Write  out  input  RHS 

computer  RHS  using  solution 
difference  of  RHS'S 

absolute  value  of  relative  differences  of  RHS 
average  absolute  relative  difference 


12.  n  —  n  +  1 

M,  —  (M  ,)gl/N 
l,n  '  l,n-l/b 

13.  Return  to  step  5  until  n  exceeds  N 

14.  Write  out  F.  matrix 

J.n 

15.  CaU  PLOT  routine 
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PLOT  ROUTINE 

1.  Input  from  LP  SOLVER:  HEADING,  N,  N  ,  F,  and  M  matrices 

2.  Input  from  IN  1  (if  used):  Ah,  H,  N  f,  andM  arrays 

3.  AM=j- 

4.  Write  out  heading 

5.  Label  vertical  axis  F(M) 

6.  Label  horizontal  axis  M 

f 

7.  versus  for  n  =  1, - N  ^ 

8.  Plot  F.  versus  M.  for  n  =  1, — N  and  i  =  1, — ,N.. 

J,n  j,n  *  J  M 

9.  STOP 
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SECTION  V 
CONCLUSION 

To  date,  the  linear  programming  method  seems  to  be  one  of  the  most 
promising  schemes  for  obtaining  the  molecular  weight  distribution  from 
sedimentation-diffusion  equilibrium  data.  There  are  five  main  features  of 
this  investigation  which  are  worthy  of  mention. 

1.  The  linear  programming  method  has  been  found  to  give  acceptable 
results  for  the  case  of  experimental  data  obtained  at  one  angular  velocity. 

2.  It  was  found  that  if  a  normal  random  error  of  the  experimental 
gradient  curve  was  about  20%  the  linear  programming  method  produced  a  MWD 
with  satisfactory  precision.  However,  if  the  experimental  error  was  weighted, 
i.  e. ,  the  concentration  gradient  curve  was  distorted  from  the  true  curve,  a  1 
or  2%  error  led  to  an  absurd  molecular  weight  distribution.  This  agrees  with 
the  findings  of  Lee  (Reference  5),  Tikhonov  and  Glasko  (Reference  12),  and 
Tikhonov  (References  13  and  14). 

3.  This  investigation  did  not  involve  any  modification  of  the  LP  solver 
routine.  The  LP  solver  limitations  were  manifested  by  spurious  points  some¬ 
times  appearing  in  the  determined  molecular  weight  distribution.  In  general, 
the  linear  programming  method  presents  only  an  overall  molecular  weight 
structure.  Therefore,  when  twenty  sets  of  molecular  weights  were  sampled, 
it  was  obvious  when  one  point  was  completely  illogical. 

4.  A  great  improvement  was  achieved  by  solving  for  F(M)  directly 
(Equation  20),  rather  than  via  f  (Equation  9).  In  the  former  case  the  matrix 
presented  to  the  LP  solver  routine  was  not  so  ill-conditioned  (Tables  I  and  II). 

5.  The  following  test  was  made  after  each  call  of  the  LP  solver  routine. 
By  having  the  calculated  MWD,  the  computer  program  could  calculate  new 

U-values  (u  cajc  )«  The  difference  between  and  Uc^c  was  then  printed. 

Also,  the  absolute  relative  differences,  the  averaged  absolute  relative  error 
for  one  set  of  molecular  weights,  and  the  averaged  absolute  relative  error 
averaged  over  all  molecular  weight  sets  were  printed.  In  general,  a  good 
"fit"  between  the  assumed  MWD  and  the  derived  MWD  showed  low  values  for 
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all  the  above  error  analyses;  however,  the  converse  was  not  always  found  to 
hold  true.  At  present  much  effort  is  being  focused  on  determination  of  a  pro¬ 
cedure  for  obtaining  a  one-to-one  correspondence  between  the  error  analysis 
criteria  and  the  "fit"  of  the  derived  MWD. 
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TABLE  I 

EXPONENTS  OF  ELEMENTS  IN  THE  MATRIX  PRESENTED  TO 
LP  SOLVER  ROUTINE  USING  EQUATION  9* 


-2 

-  | 

-  | 

- 1  0 

0 

0 

0 

0 

0 

-2 

-I 

-I 

-1  0 

0 

0 

0 

0 

0 

-2 

- 1 

-  | 

-1  0 

0 

0 

0 

0 

0 

-2 

-| 

- 1 

-1  C 

0 

0 

0 

0 

0 

-2 

- 1 

-| 

-1  0 

0 

0 

0 

0 

1 

-1 

-l 

0 

0  0 

0 

0 

0 

-2 

-5 

-1 

-  | 

0 

0  C 

0 

0 

0 

-1 

-3 

-1 

- 1 

0 

0  C 

0 

0 

0 

0 

-1 

-1 

-| 

0 

0  0 

0 

0 

0 

0 

0 

-1 

- 1 

0 

0  C 

0 

0 

1 

1 

2 

*  Value  of  determinant  =  - 

0.3932972  E -31 

TABLE  H 

EXPONENTS  OF 

ELEMENTS 

IN  THE 

MATRIX 

PRESENTED 

TO 

LP 

SOLVER  ROUTINE  USING  EQUATION 
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1 

1 

2 

3 

3  4 

4 

5 

5 

5 

1 

2 

3 

3  4 

4 

5 

5 

5 

1 

1 

2 

3 

3  4 

4 

5 

6 

6 

1 

1 

2 

3 

3  4 

5 

5 

6 

6 

1 

1 

2 

3 

3  4 

5 

5 

6 

7 

1 

2 

3 

3 

4  4 

4 

4 

3 

0 

1 

2 

3 

3 

4  4 

5 

5 

4 

2 

1 

2 

3 

3 

4  4 

5 

5 

5 

4 

1 

2 

3 

3 

4  5 

5 

6 

6 

6 

1 

2 

3 

3 

4  5 

5 

6 

7 

8 

*  Value  of  determinant  =  -0.1500715  E  17 
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EXFC.001 

r 

EXEC. 002 

r 

MAIN  PROGRAM  fxfcutivf  control 

FXFC.003 

r 

INITIALIZF  PLOTTING  ROUTINFS 

EXEC. 004 

r 

INITIALIZE  PLOT 

COUNT 

EXEC. 005 

r 

DETERMINE  SEQUENCE  OF  SUBROUTINES  CALLED 

EXEC. 006 

C 

TERMINATE  PLOTTING  BEFORE  EXITING 

EXEC. 007 

C 

WRITE  HEADING  AND  PLOT  COUNT 

EXEC. 008 

r 

EXEC. 009 

COMMON  /PLTR/  PDATAI438),  IPLTS,  HEAD (12) 

EXEC. 010 

COMMON  /PRTCTL/  WRITE 

EXEC. Oil 

r 

EXEC, 012 

LOGICAL  WRITE 

EXEC. 013 

r 

EXEC. 014 

IPLTS  =  0 

FXEC.015 

CALL  PLOTS  t  PDATA, 

438  ) 

EXEC. 016 

i 

RFAD  (5*500)  HEAD 

EXEC. 017 

WRITE  (6*600)  HEAD 

EXEC. 018 

CALL  INI  (  $900  ) 

EXEC. 019 

IF  (  WRITE  )  WRITE 

(6,600)  HEAD 

EXEC. 020 

CALL  LPS  (  $900  ) 

EXEC. 021 

CALL  PLOTR 

EXEC. 022 

GO  TO  1 

EXEC. 023 

900 

WRITE  (6,601)  IPLTS 

EXEC. 024 

CALL  PLOTE 

EXEC. 025 

CTOP 

FXFC.n?6 

600 

FORMAT  t  1246  ) 

FXEC.027 

600 

FORMA  T  (  1H1  *  12A6 

/  (  IX,  1 2  A6  )  ) 

EXEC. 028 

601 

FORMAT  (  1H0 ,  12.  17H  PLOT(S)  COMPLETE  ) 

EXEC. 029 

FNO 

EXEC. 030 

26 
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SIBFTC  RIN1.  DECK  BIN1. 

C  B  I N 1 1 

BLOCK  DATA  SUBPROGRAM  TO  SUPPLY  INPUT  DATA  WHEN  INPl.  IS  BY-PASSEDB INI. 

3  I N 1  ■ 

aLOCK  DAT  A  B  I N 1 , 

BIN1. 

COMMON  /BINT/  NW,  NX,  XL<20),  Z(20),  U(20,2D)  BIN1. 


DATA  NW  / 5 / 

DATA  NX  / 5 / 

DATA  XL  /  2.5E-6,  10.E-6,  40.E-6,  160. E-6,  640. E-6,  15*0.  / 


/ 1  .0, 


.187,  .405,  .37 2,  .144,  .005, 

.208,  .554,  .635,  .337,  .026, 

.232,  .799,1.346,  .959,  .164, 

.260,1.237,  0.  ,  0.  ,  o.  . 

.2 94,2.121,  o.,  0.,  o.. 


15*0.  / 

15*0. , 
15*0. , 
15*0. , 
15*0. , 
15*0.  ■ 


•2  00*0.  / 


B  I  N 1 . 0 1 9 
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'  BOR. 

BDR.0001 

BLOCK 

DATA 

BDR.0002 

X  IS 

R**2 

BDR.0003 

BDR.0004 

COMMON  /BLOCKS/  X ( 20  > 

BDR.0005 

X  /  36.548312, 

BDR.0006 

39.358308  , 

RDR.0007 

42.272363  , 

RDR.nnoa 

46.838556, 

BDR.0009 

48.417651  , 

BDR.0O1O 

50,012759  , 

BDR.0011 

14*0.  / 

BDR.0012 

END 

BDR.0013 
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■STOFTf  RO.  HFCK 

r 

r  RLOrK  DATA  61  ipppor-RAM  FOR  122  CO  VALUES  AND  7  W 2  V  ALUF c 

r 

BLOCK  DATA 
r 

COMMON  /RLOCKW/  W2I20) 

COMMON  /BLOCKC/  CO ( 160) 

r 

DATA  W2 

1/416930.0,  583702.0,  817182.8,  1144086.,  1601678.,  2242349., 
O-ji^opoo.,  V*#n.  / 


0  A T  A  rn 


1  / 

. 170F-5, 

. 350E-5 , 

.  530F-4 , 

.1 15E-3, 

•  850E-4 , 

.5505-4, 

.  3 1 OE-4 , 

2 

•  750E-5  , 

. 700  E-5 , 

•  650E-5 , 

.8205-5, 

. 100E-4 , 

.  135E-4, 

.  1 70E-4 , 

3 

• 260E-4 , 

. 340E-4, 

. 360E-4 , 

• 370E-4 , 

• 300E-4 , 

•  2  2  5E-4 , 

.21 5  E-4 , 

4 

.205E-4, 

.301E-4, 

. 400E-4 , 

• 435E-4, 

•  4705—4 , 

.4105-4, 

• 350E-4 , 

5 

• 280E-4 , 

.21 OE-4 , 

.  2 2 5 E-4  , 

. 23  EE-4 , 

. 2455-4 , 

.  2  60  E- 4 , 

.25  8  E-4 , 

6 

.255E-4, 

. 24  5  E— 4 , 

. 2  35E-4 , 

. 2255-4, 

.2105—4, 

.2055-4, 

• 200E-4 , 

7 

.31 OE-4 , 

•  420  5-4 , 

« 440  E— 4 , 

.4205-4, 

.3805-4, 

.  340 E-4 , 

. 300E-4 , 

8 

• 260E-4 , 

.255  E-4  » 

• 250E-4 , 

• 245E-4 , 

.2405-4, 

.232 E-4, 

.  22  5E-4 , 

9 

.21 5E-4, 

.2  05  E-4 , 

•190E-4, 

.1705-4, 

.  1635—4 , 

.1555-4, 

.17  5  E-4 , 

X 

. 1 90F-4, 

.2265-4, 

•260R-4 , 

•  24^5-4 , 

.2205-4, 

.1955-4, 

•  1 9  OE-4 , 

1 

. 1 825-4, 

.  17SF-4, 

.1706-4, 

. 1 666-4, 

•  1 60F-4 , 

.1505-4, 

.1405-4, 

2 

•  1 30F-4 , 

.  1 1 OE-4, 

.  900E-5 , 

. 650E-5 , 

.4505-5, 

.5005-5, 

•  550E-5 , 

3 

. 820E-5  , 

.1105-4, 

.1135-4, 

.1 15E-4, 

.9205-5 , 

. 700  F- 5 , 

.  60  0E- 5  » 

4 

.50OE-5, 

.400E-5, 

.3005-5, 

.2305-5, 

.  150E-5, 

.  1005-5, 

.  150E-5, 

5 

• 22°E-5 , 

.  3O0E-5 , 

.  400E-5 , 

•  500E-5  , 

.  470E-5 , 

.  4  5  0  E- 5  » 

.  420E-5 , 

6 

• 400E- 5 , 

•  5O0E-5 , 

. 600E-5 , 

. 680E-5 , 

.  750E-5 , 

. 700  E- 5  , 

•  650E-5 , 

7 

• 580E-5 , 

.  500E-5 , 

. 420E-5  , 

• 350E-5 , 

•  300E-5  * 

. 2  00E-5  » 

.  150E-5 , 

8 

.10OE-5, 

.5005-6, 

.  0  ,  0  .  ,  0 

.,  36*o . 

FND 


BD. 00001 
BD. 00002 
BD. 00003 
BD. 00004 
BD. 00005 
BD. 00006 
BD. 00007 
BD. 00008 
BD. 00009 
BD. 00010 
BD. 00011 
BD. 00012 
BD. 00013 
BD. 00014 
BD. 00015 
BD. 00016 
BD. 00017 
BD. 00018 
BD. 00019 
BD.00C2C 
BD. 00021 
BD. 00022 
BD. 00023 
BD. 00024 
BD. 00025 
BD. 00026 
5D. 00027 
BD. 00028 
BD. 00029 
BD.00C3C 
BD. 00031 
3D. 00032 
BD. 00033 
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c t bcTC  RD.  nrrx 

r 

r  nLor<  DATA  SUBPROGRAM  FOR  160  CO  VALUES  AND  7  W2  VALUES 

r 

BLOCK  DATA 
r 

COMMON  /BLOCKW/  W2 ( 2 0 ) 

COMMON  /BLOCKC /  CO (160) 
r 

DATA  W2 

1/416930.0,  583702.0,  817162.8,  1144056.*  1601678.,  ? 242349., 
7S13Q28R.,  17*0.  / 


RATA  (''R(I  ), 1  =  1,114) 
l/,6BRnRF-6 , .onnonF-E  , 
2  «23250E-4».28000F-4, 
8  .67250F— 4,.7R600E-4, 

4  ,17700E-3*.20600E-3, 

5  .42525E-3, .48775E-3, 


11000r-4 ,.lR?5nE— 4*.16O00E-4,.l 9?  50E-4 , 
33500E-4 , ,40000E-4,.4  7750E-4,  .567  50E-4, 
93750E-4, .1 102  5  E-3 , .1 2925E-3».15150E-3» 
23925E-3  , . 27 72 5E-3 » . 3 20 50E-3  ,  . 36975E-3 » 
55825E-3  ,. 6372 5E-3 ,. 72  5 75C-3  ,. 82450C-3  , 


BD. 00001 
BD. 00002 
BD. 00003 
BD. 00004 
BD. 00005 
BD. 00006 
BD. 00007 
BD. 00008 
BD. 00009 
BD. 00010 
BD. 00011 
BD. 00012 
BD. 00013 
BD. 00014 
BD. 0001 5 
BD. 00016 
BD. 00017 
BD. 00018 
BD. 00019 


6  . 9342  5E-3 , 

7  .  1 8755E-2  » 

8  .34420E-2, 
9.477325F-?  , 
X.R85R5RE-2, 

1  •  o  1 240075  , 

2  •  o  1  588050  , 
3.01 858775  , 
4.01988475  , 
5  •  R 194427  *  , 
6.01737525  , 
7.01419175  * 
8.01059475  , 
9.R0722900  , 

RATA  (CO(I) 
1/. 00450825, 

2. R02569500, 

3 .  °0 1 338  7 50  , 
.6372  5RF— 3 , 
. ?  7  7?  50F  — ^ , 

6.1  1 0  2  5 Rp- 5 , 
7 • 400000E-4 , 
8.1325R0E-4, 


.  105575F-2 

.  208825E  — 2 

•377525E-2 

•623825p-2 

. 942 1  OOF  — 2 

.01 300375 

. R 1 64047  5 

.01691575 

.0)008460 

.01920150 

• R 1 690400 

.01360175 

.01000325 

.00672375 

,1=115,160 

.00413075 

.00231925 

.0011 9025 

.658250F-S 

,?79?50F-7 
.R3750RF-4 
.  335000F-4 
.  1  1  000RE-4 


•  1 19P25F-2 
.  231925E-2 
•413R75E-2 
.672375E-? 
.01000325 
.01360175 
.0169R400 
.01920190 
.01993450 
.01891575 
.01640475 
.01300375 
.00942010 
.00623825 

.00377525 

.00208825 

.00105575 

•  487750E-?. 
.  2060008-3 
.7O6OR0E-4 
.  2  80000  5-4 
. 900000E-5 


• 133875E-2 

• 256950E-2 

. 450825E-2 

. 72  2900E-? 

.OI059475 

.01419)75 

.01737525 

.01944275 

.01088475 

.01858775 

.01588050 

.01240075 

.00885050 

.00577325 

.00344200 

.00187550 

•  9342  50E-3 
•425250F-3 
. 177000E-3 
. 672  500E-4 

•  2  32  500E-4 
• 65OQ00E-5 


.15C175F-2 

• 283950E-2 

.490  7  7  5E-2 

•  7  75275F-2 

.01119325 

.01477075 

.01781500 

.01963800 

.01978575 

.01822000 

.01533475 

.01179625 

.00829375 

.00532950 

.00313025 
.00168025 
. 82 4 500 E- 3 
. 369750F-3 
.151 500E-3 
.  567  5  OOF-4 
.  1 92  500E-4 
/ 


. 168025E-2 
•313025E-2 
• 532950E-2 
. 8  293  75  E— 2 
.01179625 
.01533475 
.01822000 
.01978575 
.01963800 
.01781500 
.01477075 
.01119325 
.007-75275 
.00490775 

.00283950 
.00150175 
• 725750E- 
.320050E- 
•  1 292  50E- 
.477500E- 
. 160000E— 


,BD. 00020 
,30.00021 
,BD. 00022 
,BD. 00023 
,BD. 00024 
,3D. 00025 
,BD. 00026 
,BD. 00027 
,BD. 00028 
,BD. 00029 
,30.00030 
,BD. 00031 
,BD. 00032 
/BD. 00033 
BD. 00034 
,BD. 00035 
,BD. 00036 
,BD. 00037 
,BD. 00038 
,BD. 00039 
,BD. 00040 
,BD. 00041 
BD. 00042 


F  NO 


BD.C0043 
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STBFTC  BD.  DECK 

r 

r  BLOCK  DATA  SUBPROGRAM  FOR  32  CO  VALUES  AND  7  W2  VALUES 


r 


r 


C 


RLOCK  DATA 

COMMON  /RLOCKW/  W2C20) 

COMMON  /BLOCKC/  C0(160) 

DATA  W2 

1/416930.0*  583702.0,  817182.8,  1144056.,  1601678.,  2242349., 
23139289.,  13*0.  / 


DATA  CO 

1 / • 91 324nOE-3 , 

2.. 2051750F-1 , 

3  ,  • 4  38 3562 E-l , 

4  ,  . 5382040E-1 , 

5 . .  5144596E  — 1  , 
6  ,  .  2496 1 94F-1 , 

7 . .  432  2  599E-2 , 
FND 


3044140F -2 , 
7538812E-1 , 
4706240E-1 , 
5461 1P8E-1  , 
4858  44  8E-1  , 
196D426E-1  , 
8 140020E-4 , 


. 6392700F-2 , 
.3041 964F-1 , 
. 4968036E-1 , 
.5467276E-1 , 
.440791  4E-1 , 
•  1 503806E-1 , 
128*0.  / 


1  o  5  3  2  72  E-l , 
3494672E-1, 
5 1 44  59  6E- 1  , 
■=>41  PF70E-1  , 
3835616E-1  , 
11 20244E-1  , 


14977 1 6F- 1 
3975646F-1 
5  2968  04E- 1 
S334344F-1 
3117200E-1 
74277,'0E-2 


BD 

BD 

BD 

BD 

BD 

3D 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 

BD 


OCOOl 

00002 

00003 

00004 

00005 

00006 

00007 

00008 

00009 

00010 

00011 

00012 

00013 

00014 

00015 

00015 

00017 

00018 

00019 

00020 

00021 

00022 
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MBMAP 

RAND  • 

100, DECK 

RAND. 001 

# 

GENERATES  UNIFORM  RANOOM  numbfps 

RAND. 002 

R=FLRAN(Y),  Y  DUMMY 

GIVES  RFAL  NUMBER 

RAND. 003 

-* 

CALL  SAVF(Z)  G1VFS 

LAST  OCTAL  VALUE 

RAND, 004 

* 

CALL  VALUE (7)  GIVES 

STARTING  OCTAL  VALUE 

RAND. 005 

FNTRY 

FLRAN 

BAND. 006 

FNTRY 

SAVE 

RAND. 007 

FNTRY 

VALUE 

RAND. 008 

FLRAN 

LDQ 

RANDOM 

RAND. 009 

MPY 

GENERA 

RAND. 010 

STO 

random 

RAND. Oil 

CL  A 

AAA 

RAND. 012 

LCL 

78 

RAND. 013 

FAn 

AAA 

RAND. 014 

TEA 

1  ,4 

RAND. 015 

V&Lf  IF 

TLA* 

8,4 

RAND. 016 

STO 

RANDOM 

RAND. 017 

TRA 

1  *4 

RAND. 01 8 

SAVF 

CLA 

RANDOM 

RAND. 019 

STO* 

3,4 

RAND. 020 

TRA 

1,4 

RAND. 021 

PANDOM 

OCT 

343777244615 

RAND. 022 

AAA 

OCT 

172000000100 

RAND. 023 

rcKjrp  A 

OCT 

34377724461  6 

RAND. 024 

CMA 

RAND. 025 
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trnFir  rlotr.  opch 

SUBROUTIMF  PLOTR 
COMMON  SCH { 4 ) *  SCF ( 4 ) 

COMMON  /PLTR/  PDATA(438),  IPLTS,  HEAD( 12 ) 

COMMON  /B I NP 1 /  DH,  H,  NC,  F<162),  CM(162)  ,  SKIP2 
COMMON  /BLPS/  N,  NM  ,  BF(IOOO),  BM(lCOO),  K. 

COMMON  /BIN2/  USEIN2 

DATA  HTITLE  / 1HM/ ,  FTITLE  '4HF(M)/ 

LOGICAL  SKIP2 
LOGICAL  USFTN2 
T  R  (  SKIP?  )  RETURN 

(  .NOT.  LKEIN2  )  GO  TO  7 
RCH ( 1  I  =  RM ( 1 ) 

SCH ( 2 )  =  BM ( K ) 
on  to  R 

7  DM  =  DH/H 

nn  5  T  =  1  ,  NC 
5  F ( I )  =  F  (  I  )  /  DM 

SCH  ( 1  )  =  AMINK  CM  (  1  )  ,  BM  ( 1 )  ) 

SCH ( 2 )  =  AMAX1 (  CM(NC),  BM ( K I  ) 

8  SCF  1 1  )  =  0. 

SCF ( 2 )  =  0. 

IP  (  USFIN2  )  GO  TO  15 

"0  1 A  T  =  l,MO 

I"  RCFI2I  =  AMAX1 (  SCF  (  2  )  ,  F (  I  >  ) 
ic  no  ?n  r  =  1  ,Y 

20  SCF ( ? )  =  AMAX1  (  SCF (  2  )  »  BF ( I  1  ) 

CALL  SCALF  (SCH,  15.,  2,  1,  10.  ) 

CALL  SCALF  (SCF,  10.,  2,  1,  10.  ) 

IF  (  USE  I N2  )  GO  TO  25 
CM ( NC+1 )  =  SCH ( 3 ) 

CM ( Nr+2 )  =  SrH<4) 

F(NC+1)  =  SOP(S) 

F ( NO+? )  =  S  C  F ( 4 ) 

25  BM(<+1)  =  SCH ( 3 ) 

RM(K+2)  =  SCH ( 4 ) 

°f(K+I)  =  sor  {  3  I 
Rp ( K+2 )  =  SCP ( 4 ) 

CALL  PLOT  (  S. ,  -11.,  — R  ) 

CALL  PLOT  (0.,  .5,-3) 

CALL  AXIS  (  0.,  0.,  HTITLE,  -1,  16.,  0.,  SCH(3),  SCH<4),  10.  ) 
CALL  AXIS  (  0.,  0.,  FTITLF,  4,  10. ,90.,  SCF(3)»  SCF ( 4 )  ,  10.  ) 
CALL  SYMBOL  (  1.,  9.5,  .25,  HEAD,  0.,  72  ) 

IF  (  .NOT.  USE  I N2  )  CALL  LINE  (  CM,  F,  NC,  1,  0,  0  ) 

CALL  LINE  (  RM,  BF,  K,  1,  1,  1  ) 

CALL  PLOT  (  15. ,  0.  ,  -3  ) 

IPLTS  =  IPLTS  +  1 

RFTURN 

CMD 


PLOTR.Ol 

PLOTR.02 

PLOTR.03 

PLOTR.04 

PLOTR.05 

PLOTR.06 

PLOTR.07 

PLOTR.08 

PLOTR.09 

PLOTR.IO 

PLOTR.ll 

PLOTR.12 

PLOTR.13 

PL0TR.14 

PLOTR.15 

PLOTR.16 

PLOTR.17 

PLOTR.18 

PLOTR.19 

PLOTR.20 

PLOTR.21 

PLOTR.22 

PLOTR.23 

PLOTR.24 

PLOTR.25 

PLOTR.26 

PLOTR.27 

PLOTR.28 

PLOTR.29 

PLOTR.30 

PLOTR.31 

PLOTR.32 

PLOTP.33 

PLOTP.34 

PLOTR.35 

PLOTR.36 

PLOTR.37 

PLOTR.38 

PLOTR.39 

PLOTR.40 

PLOTR.41 

PLOTR.42 

PLOTR.43 

PL0TR.44 

PLOTR.45 

PLOTR.46 

PLOTR.47 

PLOTR.48 

PLOTR.49 
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■MPFTr  ORDER.  DECK 

SUBROUTINE  ORDER  (  X,  Y.  NM,  N *  K  ) 

COMMON  /PRTCTL /  WRITE 
DIMENSION  X(l) *  Y ( 1  ) 

LOGICAL  WRITE 
J  =  0 

DO  ] n  LN  =  1 *N 
M?0  =  (LN-1)  *  20 

1  =  +  i 

J  =  J  +  1 

X(J)  =  X  C  L ) 

YfJ)  =  Y  C  L  1 
DO  in  I  =  2 » NM 
L  =  N20  +  I 

J  =  J  +  1 

X  (  J  1  =  X  (  L  ) 

in  v ( j )  =  YIL) 

Y  =  J 

20  TOST  =  n. 

nn  nn  t  =  2,J 

TF  (  X(I-l)  ,LF.  X(I)  I  GO  TO  SO 

XP  =  Xf  T-n 
yp  =  y(  T-n 
XU-1)  =  XU) 

YU-1  )  =  YU  ) 

XU)  =  xs 
Y(I  )  =  YS 
TUST  =1. 

30  CONTINUE 

IF  (  (TEST  .FQ.  0.)  .OR.  (  J  .ECU  2  )  )  GO  TO  40 
J  =  J  -  1 
GO  TO  20 

40  IF  (  WRITE  ).  WR I TF  (6*600)  (I,  X(I),  Y(I),  1  =  1, K') 

ncTI  |PM 

600  FORMAT  (  1 H 1 ,  15X,  1HM,  IPX,  1HF  /  1H  / 

X  (  IX,  13,  2E20.7 

END 


ORDER. 01 
ORDER. 02 
ORDER. 03 
ORDER. 04 
ORDER. 05 
ORDER. 06 
ORDER, 07 
ORDER. 08 
ORDER. 09 
ORDER. 10 
ORDFR.ll 
ORDER. 12 
ORDER. 13 
ORDER, 14 
ORDER. 15 
ORDER. 16 
ORDER, 17 
ORDER. 18 
ORDER. 19 
ORDFR.20 
ORDER. 21 
ORDER. 22 
ORDER. 23 
ORDER, 24 
ORDER. 25 
ORDER. 26 
ORDER. 27 
ORDER. 28 
ORDER. 29 
ORDER. 30 
ORDER. 31 
ORDER. 32 
ORDER. 33 
ORnFR.34 
ORDER. 35 
ORDER. 36 
ORDER. 37 
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$  I BFTC  DETA.  DECK 

SUBROUTINE  DETA  (  A,  B*  NM*  NROW  ) 
COMMON  /PRTCTL/  WRITE 
DIMENSION  A (51 *120),  R(NM,NM) 
LOGICAL  WR  T  T^ 

IF  (  WRITE  )  WRITP  (6*600) 

J  =  1 

ID  DO  20  I  =  1,NM 
DO  20  K  =  1  * NM 
J I  =  J  +  I 
20  R  (  I  *  K  )  =  A ( J I  *  K  ) 

D  =  DET  (  B.  NM  ) 

IF  (  WRITE  )  WRITE  (6*601)  J,  D 
J  =  J  +  NM 

IF  (  ( J+NM-1 I  .GT.  NROW  )  RETURN 
GO  TO  in 

600  FORMAT  (  1 3H0DFTFRM I N ANTS  I 

601  FORMAT  (  3 X  *  14*  E19.7  I 
FND 


DETA. 001 
DETA. 002 
DETA. 003 
DETA. 004 
DETA. 005 
DFTA.006 
DETA. 007 
DETA .008 
DETA. 009 
DETA. 010 
DETA. Oil 
DETA. 012 
DETA. 013 
DETA. 014 
DETA. 015 
DETA. 016 
DETA. 017 
DETA. 018 
DFTA.019 
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STBFTC  DLETE.  DECK 

SUBROUTINE  DLETE  (  IRHSi  NROW,  NRP1,  NCOL,  NM ,  *  ) 

COMMON  JRHS(IOO) 

COMMON  /BIN1/  NW,  NX,  XL(21)>  Z(20),  U(20,20) 

COMMON  /DLT/  NWW,  NXX 
COMMON  /PRTCTL/  WRITE 
LOGICAL  WRITF 
XPHS  =  ?  *  IRUR 

I F  (  KRHS  .GT.  100  )  GO  TO  900 
READ  (5*500)  < JRHS( I > ,1*1 ,KRHS) 

IF  (  WRITE  )  WRITE  (6*601)  IRHS*  ( JRHS ( I) *  I  =  1  * KRHS ) 

DO  10  I  =  1  *  KRHS  *  2 
K  =  JRHS ( I ) 

L  =  JRHS ( I +1 ) 

10  U(K,L)  =  0. 

NWW  =  NW 
!  =  0 

on  SO  K  =  1  *  NW 
USUM  =  0, 

DO  ?n  L  =  1  *  NX 
20  USUM  =  USUM  +  U(K*L) 

IF  {  USUM  .EQ.  0.  )  GO  TO  28 
1  =  1  +  1 
DO  25  J  =  1  *  NX 
25  U ( I  *  J )  =  U  (  K  *  J  ) 

XL ( I )  =  XL  (K) 

GO  TO  30 
28  NWW  =  NWW-  1 
30  CONTINUE 
NXX  =  NX 
J  =  0 

on  so  L  =  1 » NX 
USUM  =  0* 
on  40  K  =  1  *  NWW 
40  USUM  =  USUM  +  U  (  K  »  L  ) 

IF  (  USUM  .FO.  0.  )  GO  TO  48 
J  =  J  +  1 
DO  45  I  =  1  *  NWW 
45  U ( I , J  )  =  U ( 1  *  L ) 

Z(J)  =  Z(L) 

GO  TO  50 
48  NXX  =  NXX  -  1 
5 o  CONTINUE 

NROW  =  NXX  *  NWW 
NRP1  =  NROW  +  1 
NCOL  =  NM  +  ( 2  *  N  R  OW ) 

IF  (  .NOT.  WRITS  )  RETURN 
WRITE  (6,602)  NWW,  NXX 
DO  60  I  =  1 , NWW 

60  WRITE  (6,603)  I,  ( U ( I , J ) , J=1 ,NXX ) 

RETURN 

oon  WR I Tr  (6,600)  IRHS 
RETURN  1 

500  FORMAT  (  27 (  211, IX)  ) 

600  FORMAT  (  45HCNUMBER  OF  U  MATRIX  DELETIONS  GREATER  THAN  50  / 

X  7H0IRHS  =  14  ) 

601  FORMAT  (  1HA ,  13,  41H  ELEMENTS  OF  MATRIX  -U-  HAVE  BEEN  DELETED  / 

X  33 ( 2X , 2  1 1 )  ) 

602  FORMAT  (  20HAU  MATRIX  (ADJUSTED),  110,  5H  ROWS,  16,  8H  COLUMNS  ) 

603  FORMAT  (  4H0R0W,  14*  IX,  6E20.7  /  (9X,  6E20.7)  ) 

FNO 


DLETE001 
DLETE002 
DLETE003 
DLETE004 
DLETE005 
DLETE006 
DLETE007 
DLFTF008 
DLFTE009 
DLETE010 
DLETE011 
DLETE012 
DLETE013 
DLETE014 
DLETE015 
DLETE016 
DLETE01 7 
DLETF018 
DLETE019 
DLETF020 
DLETE021 
DLETE022 
DLETE023 
DLETE024 
DLETE025 
DLETE026 
DLETE027 
DLETE028 
DLETE029 
DLETE030 
DLETE031 
DLETE032 
DLETE033 
DLETF034 
DLETE035 
DLETE036 
DLETE037 
DLETE038 
DLETE039 
DLETE040 
DLETE041 
DLETE042 
DLETF043 
DLETE044 
DLETE045 
DLETE046 
DLETE047 
DLETE048 
DLETE049 
DLETEO&O 
DLETE051 
DLFTF052 
DLETE053 
DLETE054 
DLETE055 
DLETE056 
DLETE057 
DLETE058 
DLETE059 
DLETE060 
DLETE061 
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$  T  BFT C  LPS.  DECK 

SUBROUTINE  LPS  (  *  ) 

COMMON  I E ( 20  )  *  E ( 5 1  *  5 1 ) »  X<51),  P(51),  Y(51),  JH(51), 

X  KB (120),  KOUT ( 7 ) »  ERP<8) 

COMMON  / B I N 1 /  NW ,  NX,  XL(20),  Z(20),  U(20,20) 

COMMON  /RLPS/  N,  NM ,  BF(20,50),  BM(2O,50),  K 
COMMON  /DLT/  NWW,  NXX 
roMMON  /PRTrTL/  WRTTF 

DIMENSION  A ( 5 1 » 1 20 )  *  B<51),  I N  F I X ( 8  >  »  T0L(4)»  BA ( 400 ) 

D I MFNS I  ON  A  B  R ( 50)  ,  BS(51) 

LOGICAL  WRITE 

DATA  NMAX/50/,  NMMAX/20/*  NCOLMX/120/*  NROWMX/50/ 

NAME  LIST  /LP/  XM11,  G,  NM,  N,  IRHS»PERT 
PERT  =  0. 

TRHS  =  0 
RFAD  ( 5  ,  LP  ) 

IF  t  WRITF  )  WRITF  ( 6 , 6 n  0 /  XMJI,  G ,  NM,  N 
IF  (  i\j  ,GT.  MMAX  )  GO  TO  900 
IF  (  Mm  .GT.  NMMAX  )  GO  TO  904 
NjPOW  =  NW  *  NX 

IF  (  NROW  .GT.  NROWMX  )  GO  TO  902 
NCOL  =  NM  +  (2*NROW) 

IF  (  NCOL  .GT.  NCOLMX  )  GO  TO  901 
NR PI  =  NROW  +  1 
NWW  =  NW 
NXX  =  NX 

IF  (  I RHS  .GT.  0  ) 

X  CALL  DLETF  (  IRHS,  NROW,  NRP1,  NCOL,  NM ,  $90?  ) 

20  DO  ?0  I  =  1 , NRP 1 
DO  SO  J  =  1 , NCOL 
?0  A (  I , J  )  =  0. 

dd  OF  i  =  1 , NMMAX 

DO  7  F  J  =  MMAX 

39  BF ( I , J)  =  0. 

INFIX (1)  =  4 
I NF I  X ( 2 )  =  NCOL 
INFIX!?)  =  91 
I  NF  I X ( 4 )  =  NRP 1 
INFIX! 9)  =  2 
INFIX (6)  =  1 
INFIX (7)  =  900 
INFIX (8)  =  20 
Tol (1 )  =1 . F— 7 

TOL ( 2 )  =  I.F-7 

T OL ( 3 )  =  l.F-6 
TOL ( 4 )  =  l.F-lo 
DRM  =  0. 

B ( 1 )  =  0. 

DO  50  I  =  1 , NROW 
IN  =  ( I  — 1 )  /  NXX 
K  =  t  N  +  1 
L  =  I-  ( I N*NXX ) 

F(  1  +  1  )  =  !KK,L) 

IF  (  PFRT  .GT.  0.  )  B ( 1+1 )  =  B( 1+1 ) *( l.+2.*PFRT*( FLRAN (X1-.5) ) 
80  f S ( 1  +  1  )  =  B ( 1  +  1) 

IN  =  1 

AORTG  =  SORT!  G  ) 

GMG  =  SORT G  -  (1.  /  SQRTG) 

XN  =  1.  /  FLOAT ( N ) 

BM(  1  ,1  )  =  XMU 

NROWS  =  NROW 


LPS. 0001 
LPS. 0002 
LPS. 0003 
LPS. 0004 
LPS. 0005 
LPS. 0006 
LPS. 0007 
LPS. 0008 
LPS. 0009 
LPS. 0010 
LPS. 0011 
LPS. 0012 
LPS. 0013 
LPS. 0014 
LPS. 0015 
LPS. 0016 
LPS. 0017 
LPS. 0018 
LPS. 0019 
LPS. 0020 
LPS. 0021 
LPS. 0022 
LPS. 0023 
LPS. 0024 
LPS. 0025 
LPS. 0026 
LPS. 0027 
LPS. 0028 
LPS. 0029 
LPS. 0030 
LPS. 0031 
LPS. 0032 
LPS. 0033 
LPS. 0034 
LPS. 0035 
LPS. 0036 
LPS. 0037 
LPS. 0038 
LPS. 0039 
LPS. 0040 
LPS. 0041 
LPS. 0042 
LPS. 0043 
LPS. 0044 
LPS. 0045 
LPS. 0046 
LPS. 0047 
LPS. 0048 
LPS. 0049 
LPS. 0050 
LPS. 0051 
LPS. 0052 
LPS. 0053 
LPS. 0054 
LPS. 0055 
LPS. 0056 
LPS. 0057 
LPS. 0058 
LPS. 0059 
LPS. 0060 
LPS. 0061 
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55  IF  (  WRITE  )  WRITE  (6*612)  LN.  N 
NROW  =  NROWS 
WRP1  =  NROW  +  1 
DO  fen  J  =  1,NM 

RW(J,LN)  =  BM(1,LN)*  G**  <  J-l  ) 

DRM  =  GMG  *  RM(J»LN) 

DO  6^  I  =  1  *  NROW 
TN  =  ( I  — 1 )  /  NXX 
K  =  IN  +1 
L  =  I-  (  I N*NXX ) 

XLKM  =  XL ( K )  *  BM(J»LN) 

XNUM  =  XLKM**2  *  EXP (  -XLKM  *  Z(L)  ) 

DENOM=  1.  -  EXP  f  -XLKM  ) 

60  A  U  +  l  »  J  )  =  (XNUM  /  DENOM)  *  DBM 
IT  =  1 

DO  64  I  =  2*NRP1 
BCI!  =  BS(  I) 

IB  C  BUS  .FQ.  0.  )  00  TO  64 

11=11+1 
BCII)  =  B!  I  ) 

DO  62  J  =  l.NM 
62  AUI(J)  =  A  (  I  *  J  ) 

64  CONTINUE 
MRP1  =  II 
NROW  =  NRP 1  -  1 
INFIX (4 )  =  NRP 1 
J  =  NM 

DO  68  I  =  2  » NRP 1 
J  =  J  +  1 
A(ltJ)  =  I .0 
All *J)  =  1.0 
J  =  J  +  1 
A ( 1 ♦ J )  =  1.0 
68  A ( I  * J  )  =  -1.0 

CALL  DETA  (  A,  BA,  NM,  NROW  ) 

IF  (  WRITE  )  WRITE  (6,601)  (M,M=1,30) 

DO  80  I  =  2,NRP1 
DO  70  J  =  1 , NM 
T  r ! J )  =  -9  9 

70  IE  (  A ( I , J )  .NF.  0.  )  I F ( J )  =  AL0G1 O (  A ( I , J )  ) 

TM1  =  T  -  1 

80  IE  (  WRITF  )  WRITE  (6,602)  IM1,  (I F ( J ) , J= 1 , NM ) 

CALL  SIMPLX  (  INFIX,  A,  B,  TOL ,  PRM ,  KOUT,  ERR,  JH,  X,  P,  Y ,KB*E 
IF  (  WRITE  )  WRITE  (6,6C3) 

DO  90  J  =  1 , NRP 1 
MX  =  JH ( J ) 

90  IF  (  (MX  • GT .  0)  .AND.  (MX  .LE.  NM )  )  BF { MX  *  LN )  =  X(J! 

ARR(LN)  =  0. 

DO  110  I  =  2 , NRP 1 
sr  =  (i, 

DO  100  J  =  1,NM 

]Of>  BC  =  RC  +  RE  (  J  ,  L  N !  *  A  (  T  ,  J  ) 
so  =  p ( t )  -  pr 
AB  =  ABS (  BO  )  /  B< I ) 

TE  (  WRITE  )  WRITD  (6,604)  BID,  BC,  90*  AB 
110  ARR(LN)  =  ARR ( LN )  +  AB 

ABR(LN)  =  ABR(LN)  /  FLOAT (  NROW  ) 

IF  (  .NOT.  WR I TF  )  GO  TO  115 
WRITE  (6,605 )  ABR( LN) 

WRITE  (6,606)  LN,  ( BF ( M , LN ) ,M= 1 , NM ) 

WRITE  (6,610) 


LPS.0062 
LPS.0063 
L  PS • 0064 
LPS.0065 
LPS.0066 
LPS.0067 
LPS.0068 
LPS.0069 
LPS.0070 
LPS.0071 
LPS.0072 
LPS.0073 
LPS.0074 
LPS.0075 
LPS.0076 
LPS.0077 
LPS.0078 
LP5.0079 
LPS.0080 
LPS.0081 
LPS.0082 
LPS.0083 
LPS.0084 
LPS.0085 
LP  S. 0086 
LPS.0087 
LPS.0088 
LPS.0089 
LPS.OORO 
L PS. 00 91 
LPS.0092 
LPS.0093 
LPS.0094 
LPS.0095 
LPS.0096 
LPS.0097 
LPS.0098 
LPS.0099 
L  PS  .  01 00 
LPS.0101 

LPS. 01 o? 

LPS.0103 
) LPS • 0104 
LPS.0105 
LPS. 0106 
LPS, 0107 
LPS. 0108 
LPS. 0109 
LPS. 0110 
LPS. 0] II 
LPS. 01 12 
LPS. 0113 
LPS. 01 1 4 
LPS. 0115 
LPS. 0116 
LPS. 0117 
LPS, 0118 
LPS. 0119 
LPS. 0120 
LPS, 0121 
LPS. 0122 
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113  |_  W  =  I  M  +  1 

Tc  (  l_M  ,GT.  N  1  r-n  TO 
RM(i,lN)  =  ^MdiLN-l)  *  G**XN 
on  TO  3  3 

120  IF  (  WRITE  )  WRITE  (6.610) 
WRTTF  (6.616) 

ARBAVG  =  0. 

DO  130  I  =  1 ,N 

ARBAVG  =  ARBAVG  +  ARR(T) 

130  WRITF  (fa.614)  I.  ABR(I),  (BF 
ARBAVG  =  ARBAVG  /  FLOAT (  N  ) 

WRITE  (6,613)  ARBAVG 

C  ALL  ORDER  (  BM,  BF  ,  NM  ,  N ,  K 

RFTURN 

orr  WRITF  (6.607)  N 
RETURN  1 

00]  WRITE  (6,608)  NCOL 
RFTURN  1 

902  WRITE  (6,609)  NROW 

on?  ROT! ! R N  ) 

904  WRITE  (6,611)  NM 
RFTURN  1 

600  FORMAT  (  33HAFIRST  MOLECULAR 


( BF ( M . I ) ,M= 1 , NM ) 

N  ) 


33HAFIRST  MOLECULAR  WEIGHT  =  E16.7 

33H  MOLECULAR  WEIGHT  MULTIPLIER  =  E16.7 

33H  NUMBER  OF  MOLECULAR  WT.  VALUES  =  18 

33H  NUMBER  OF  LP  SETS  FOR  SOLUTION  =  18 


601  FORMAT 
X 

607  FORMAf 

603  FORMAT 
X 

604  FORMAT 

605  FORMAT 

606  FORMAT 

607  FORMAT 

608  FORMAT 
XON  / 

609  'FORMAT 
X 

610  format 

611  FORMAT 


!  1  HA ,  57X,  1 4HA  MATRIX  (LOG)  /  1H0,  60X,  7HC0LUMNS  / 
5H  ROWS,  3014  /  1H  ) 

!  IX,  13,  IX,  3014  /  (  5  X  ,  3014  )  ) 

:  1H1 ,  7X »  9HINPUT  RHS ,  7X»  1 2HC0MPUTED  RHS,  6X» 

14HRHS  DIFFERFNCE,  4X,  16HABS  REL  D1FF  RhS  / 
1  IX,  4E18.7  ) 

[  34H0AVERAGE  RELATIVE  DIFFERENCE  RHS  =  E16.7  ) 

[  9HOSOLUT I  ON ,  14,  7X,  7E16.7  /  (20X,  7E16.7)  ) 

1  4H0N  =  14,  42H  IS  GREATER  THAN  DIMENSION  FOR  NO.  OF 
[  56HONI imrfR  OF  COLUMNS  FOP  -A-  MATRIX  GREATER  THAN  DI 
7HOMROW  =  T  4  ) 


63WONUMBFR 
/  7HOROW  = 

1  HI  ) 

5  H  0  N  M  =  14. 


OF  ROWS  FOP  -A-  MATRIX  GRFATFR  THAN  DIMEn 


GREATER  THAN  DIMENSION  FOR  NUMBER 


XOLUTTONS  PER  SET 


612  FORMAT 

613  FORMAT 
X 

614  FORMAT 

615  FORMAT 

X  / 


1H0 ,  57X »  3HSET  13,  3H  OF  13  ) 

1H0  /  20X ,  38HTHE  AVERAGE  REL.  DIFF.  FOR  ALL  SETS 
E16.7  ) 

1H0 ,  14,  E22.7,  9X ,  6F16.7  /  (36X,  6E16.7)  ) 

WO  SET ,6X , 15HAVG.  Rfl.  ERROR  ,  12X , 1 1HS0LUT IONS  •  ,f 


LPS.0123 
LPS.0124 
LPS.0125 
LPS.0126 
LPS.0127 
LPS.0128 
LPS.0129 
LPS.0130 
LPS.0131 
LPS.0132 
LP  S , 01 33 
LPS.0134 
LPS.0135 
LPS.0136 
LP  S . 0 137 
LPS.0138 
LPS.0139 
LPS.0140 
LPS.0141 
LPS.0142 
LPS.0143 
LPS.0144 
LPS.0145 
LPS.0146 
LPS.0147 
LPS.0148 
LPS.0149 
LPS.0150 
LPS.0151 
LPS.0152 
1H  ) LPS*  0 1 53 
LPS.0154 
LPS. 0155 
LPS. 0156 
SETSJLPS.0157 
MENSILPS.0158 
LPS. 0159 
SION  LPS. 0160 
LPS. 0161 
LPS. 0162 
OF  SLPS.0163 
LPS. 0164 
LPS. 0165 
LPS. 0166 
LPS. 0167 
LPS. 0168 
1H*)  LPS. 0169 
LPS. 0170 

LPS. 0171 
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*.  T  BFTO  INI.  DECK 

INI. 0001 

SUBPOUT  INF  INI  (  *  ) 

INI. 0002 

COMMON  /BLOCKS/  X(20) 

INI. 0003 

COMMON  /BLOCKC/  CO (160) 

INI, 0004 

COMMON  /BLOCK W/  W2!20) 

INI. 0005 

COMMON  /B I N 1 /  NW,  NX,  XL(20),  Z(20),  U(20,20) 

INI. 0006 

COMMON  /B I N2 /  USFIN2 

INI, 0007 

COMMON  /R I NP 1 /  OH,  H,  NC,  FI162), 

HO ( 162 ) ,  SKIP2 

INI. 0008 

COMMON  /PRTCTL/  WRITE 

INI. 0009 

r 

INI. 0010 

PAT A  NXMAX/2 0/,  NHMAX/160/,  NWMAX/20/ 

INI. 0011 

r 

INI. 0012 

LOGICAL  WRITE,  SKI  PI,  SKI ->2 

INI. 0013 

LOGICAL  USFJM2 

INI, 0014 

r 

INI. 0015 

NAMF  LIST  / I NP 1 /  NW,  NX,  HI,  OH, 

NH,  H,  CZ,  XI,  DX, 

INI .0016 

X  1ST,  LST ,  LAST,  WRITE,  SKIP1,  SKIP2,  USE  I N2 

INI. 0017 

c 

INI. 0018 

WRITE  =. FALSE. 

INI. 0019 

SKI  PI  =  .FALSE. 

INI. 0020 

f<  T  P?  =  .FALFF, 

INI .0021 

UGFIN?  =  .FALSF. 

INI  .0022 

|  AFT  =  1 

INI .002? 

pc AO  ( 5 , T  N°1 S 

INI. 0024 

IF  (  LAST  ,F0.  0  )  RFTURN  1 

INI. 0025 

IF  (  SKIP1  I  GO  TO  75 

INI. 0026 

IF  (  NX  .GT.  NXMAX  )  GO  TO  900 

INI. 0027 

IF  (  NH  .GT,  NHMAX  )  GO  TO  901 

INI. 0028 

IF  (  NW  .GT.  NWMAX  )  GO  TO  902 

INI. 0029 

TF  <  USFIN2  !  CALL  TN2  <  $85  ) 

INI. 0030 

NC  =  LST  -  1ST  +  1 

INI. 0031 

IF  (  NC  .GT.  NHMAX  )  GO  TO  903 

INI .0032 

IF  (  XI  .LE,  0.  I  GO  TO  5 

INI. 0033 

X(l)  =  X] 

INI. 0034 

DO  10  L  =  2, NX 

INI. 0035 

in 

X  C  L I  =  XI  +  FLOAT ( L  — 1 I  *  DX 

INI. 0036 

6 

HO( 1 )  =  HI 

INI .0037 

HO ! 1 )  =  HI 

INI. 0038 

nn  20  N  =  2,MH 

INI. 0039 

20 

H0(N1  =  HI  +  FLOAT ( N-l )  *  OH 

INI. 0040 

on  so  I  =  1 ,NC 

INI. 0041 

IC  =  I  +  1ST  -  1 

INI. 0042 

F  (  T  )  =  CO ( IC)  /  CZ 

INI .0043 

30 

HOE  1 i  =  HOI IC)  /  H 

INI .0044 

R  =  XI  NX) 

INI. 0045 

RMM  =  p  -  XI 

INI. 0046 

rrn  40  K  =  1,NW 

INI. 0047 

4n 

XI  (X)  =  PMM  *  m  *  w 2 ( K ) 

INI. 0048 

DO  60  L  =  1 , NX 

INI. 0049 

Ejn 

7  !  L )  =  { R  -  X ( L )  )  /  BMM 

INI. 0050 

DO  70  K  =  1 , NW 

INI. 0051 

DO  70  L  =  1 , NX 

INI. 0052 

U(K,L)  =  0. 

INI .0053 

DC  70  N  =  1 »NC 

INI. 0054 

XLM  =  XL ! K  )  *  HO { N ) 

INI, 0055 

XNUM  =  XLM**2  *  EXP (  -XLM*Z(L)  ) 

INI. 0056 

rsFNOM=  1.  -  CXP(  -XLM  ) 

■"1  U(K»L)  =  U  ( K  *  L  )  +  (  XNUM  /  DENOM  )  *  F  { N  ) 
7=  JP  (  .NOT.  WPITP  )  00  TO  85 

y9TTc  (6,601)  NW,  <W?(K) ,K=1,NW) 

WRITE  (6,60?)  NX ,  (7(U  ,!=!,  NX) 


INI, 0057 
INI. 0058 
INI .0059 
INI. 0060 
INI. 0061 
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WRITE  (6,603)  NW,  NX 
DO  80  K  =  1 » NW 

80  WRITE  (6*604)  K*  ( U ( K *L ) *L= 1 » NX ) 

85  RETURN 

900  WRITE  (6,605)  NX 
STOP 

901  WRITE  (6*606)  NH 
STOP 

902  WRITE  (6,607)  NW 
STOP 

P03  W&TTF  (6*608)  NT 
STOP 

601  FORMAT  (  1HA ,  14*  40H  VALUES  OF  ANGULAR  VELOCITY  SQUARED  -  W 2  / 

X  1H  /  (  9X ,  6E20.7  )  ) 

602  FORMAT  (  1H0*  14,  71H  VALUES  OF  THE  FUNCTION  OF  DISTANCE  SQUARED  F 

XROM  CENTER  OF  ROTATION  -  Z  /  1H  /  (  9X,  6E20.7  )  ) 

603  FORMAT  (  9H0U  MATRIX,  110,  5H  ROWS,  16,  8H  COLUMNS  ) 

604  FORMAT  (  4HOROW,  14,  IX,  6E20.7,  /  (  9X,  6E20.7  )  ) 

605  FORMAT  (  5H0NX  =  14,  33H  IS  GREATER  THAN  DIMENSION  FOR  X.  ) 

606  FORMAT  (  5H0NH  =  14,  33H  IS  GRFATFR  THAN  DIMENSION  FOR  H.  ) 

607  FORMAT  (  5H0NW  =  14,  37H  IS  GREAIER  THAN  DIMENSION  FOR  OMEGA.  ) 

608  FORMAT  (  5H0NC  =  14,  41H  IS  GREATER  THAN  DIMENSION  FOR  SFLECTFR  H) 


INI. 0062 
INI. 0063 
INI. 0064 
INI. 0065 
INI. 0066 
INI. 0067 
INI. 0068 
INI. 0069 
INI. 0070 
INI. 0071 
INI .0072 
INI. 0073 
INI. 0074 
INI, 0075 
INI. 0076 
INI. 0077 
INI. 0078 
INI. 0079 
INI. 0080 
INI. 0081 
INI. 0082 
INI. 0083 


END 


INI. 0084 
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SIBFTC  IN2.  DECK 

IN2.0001 

SUBROUTINE  IN2  (  *  ) 

IN2.0002 

r 

SCRATCH  STORSGF 

IN2.0003 

COMMON  DCDX (20*20) 

IN2.0004 

r 

IN2.0005 

DIMENSION  DNOR (20*20) 

IN2.0006 

r 

IN2.0007 

LOGICAL  WRITE 

IN2.0008 

r 

IN2.0009 

COMMON  /B  INI  /  NW,  NR*  XL ( 2  0 )  *  Z ( 20 ) *  U(20*20) 

IN2. 0010 

COMMON  /BLOCKW/  W2 ( 20 ) 

IN2.0011 

COMMON  /BLOCKR/  XI  20) 

I  M2. 00 12 

COMMON  /PRTCTLZ  WRITE 

IN2.0013 

r 

IN2.0014 

DATA  NRMAX/20/,  NWMAX/20/ 

IN2.0015 

r 

IN2.0016 

NAME  LIST  / 1 NP2/  NW,  NR,  H,  CO,  Rl,  DR*  DCDN  * '  DNDR  *  W2 

IN2.0017 

r 

1N2.0018 

RF  AD  ( 5  *  I NP2 ) 

IN2.0019 

r 

IN2.0020 

TF  (  NR  .GT ,  NRMAX  )  GO  TO  900 

IN2.0021 

IF  (  NW  .GT.  NWMAX  )  GO  TO  901 

IN2.0022 

r 

IN2.0023 

DO  10  L  =  1 , NR 

IN2.0024 

IF  (  Rl  .LE.  0.  )  GO  TO  5 

IN2.0025 

X  ( L )  =  t  Rl  +  FLOAT  (L-l)  *  DR.  )  **2 

IN2.0026 

r 

IN2.0027 

5 

DO  10  K  =  1 , NW 

IN2.0028 

10 

DCDX(K,L)  =  1.0  /  SORT ( X  t  L ) )  *  DCDN  *  DNDR ( K  *  L ) 

IN2.0029 

r 

IN2.0030 

XM  =  X ( 1 ) 

IN2.0031 

n  =  X(NR) 

IN2.O032 

r 

IN2.0033 

DO  20  K  =  1,  NW 

IN2.0034 

20 

X L  (  K  )  =  (B-XM)  *  H  *  W2  (  K  ) 

IN2.0035 

r 

IN2.0036 

DO  30  L  =  1,  NR 

IN2.0037 

30 

Z  (  L  !  =  t  B— X ( L )  )  /  (B-XM) 

IN2.0038 

r 

IN2.0039 

DO  40  L  =  1 , NR 

IN2.0040 

DO  40  K  =  1 , NW 

IN2.0041 

40 

U  t  K  ,  L  )  = ( ( B— XM )  /  CO)  *  DCDX ( K  *  L ) 

IN2.0042 

r 

1N2.0043 

IF  (  .NOT.  WR I Tc  )  GO  TO  60 

IN2.0044 

r 

I N  2 . 0045 

WRITF  (6*601)  NW,  ( W2 ( K ) ,K= 1 » NW ) 

I N  2 . 0046 

WRITE  (6*602)  NR*  (  Z ( L ) * L= 1  *  NR ) 

IN2.0047 

WRITE  (6,603)  NW,  NR 

IN2.0048 

DO  50  K  =  1 , NW 

IN2.0049 

50 

WRITE  (6,604)  K,  ( U ( K , L ) , L= 1 , NR ) 

IN2.0050 

C 

IN2.0051 

60 

RETURN  1 

IN2.0052 

900 

WRITE  (6*605)  NR,  NRMAX 

IN2.0053 

ST  OR 

IN2.0054 

00] 

WRITR  (6*606!  NW,  NWMAX 

IN2.0055 

FTOD 

I  N 2. 00  56 

r 

1N2.0O57 

601 

FORMAT  (  1HA*  14,  40H  VALUES  OF  ANGULAR  VELOCITY  SQUARED  -  W2  / 

IN2.0058 

X 

:  1H  /  (  9X,  6E20.7  )  1 

1N2.005O 

602  FORMAT  (  1H0»  14*  71H  VALUES  OF  THE  FUNCTION  OF  DISTANCE  SQUARED  FIN2.0060 
XROM  CENTER  OF  ROTATION  -l  /  1H  /  (  9X,  6E20.7  )  )  IN2.0061 
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603 

FORMAT  i 

[  9H0U  MATRIX 

,  I  10 

,  5 H  ROWS , 

16, 

8H  COLUMNS 

) 

IN2.0062 

604 

FORMAT  i 

[  4H0R0W* 

14, 

IX, 

6E20.7  /  ( 

9X» 

6E20.7  )  ) 

IN2.0063 

605 

FORMAT  i 

(  5H0NR  = 

14, 

33H 

IS  GREATER 

THAN 

DIMENSION 

FOR 

XI 

,12 

,2H)  . 

) IN2.0064 

606 

FORMAT 

FMO 

(  5H0NW  = 

14, 

33H 

IS  GREATER 

THAN 

DIMENSION 

FOR 

w< 

,12 

,2H)  . 

) IN2.0065 
IN2.0066 
I N2*  0067 
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SIBFTC  DET .  DECK  L.  B.  FALL 

FUNCTION  DET(A,N) 

C 

C  DETERMINANT  EVALUATING  FUNCTION 

C 

r  FUNCTION  DFT  <  A  *N  )  COMPUTFS  THF  DETFRM I NANT 

C  OF  THE  N-TH  ORDER  MATRIX  A,  WHICH  MUST  BE 

C  DIMENSIONED  ACN*N).  THF  ORIGINAL  MATRIX  A 

r  MOT  ALTFRFO. 

r 

r  612  CELLS  OF  BLANK  COMMON  ARF  USFD 

C 

C  TO  CHANGE  DIMENSIONS,  CHANGE  DIMENSIONS  OF  ARRAYS  B  AND  PIV, 

C  AND  ALSO  CHANGE  VALUE  OF  NMAX  IN  THE  DATA  STATEMENT. 

r 

DIMENSION  A(N,N) 

COMMON  D,  I,  II,  J,  K,  KCT ,  KFROM,  KTO,  NN,  NPR ,  RLE,  TPE 
roMMON  R(24 ,24),  P I V ( 2 A ) 

DATA  NMAX/  24/ 

r 

r  TEST  ARGUMENT  N  TO  PREVENT  OVERFLOWING  BLANK  COMMON 

NN  =  N 

IF  (  NN  .GT.  NMAX  .OR.  NN  .LE.  0  !  GO  TO  100 

MOVE  INPUT  MATRIX  A  TO  SCRATCH  MATRIX  B 

DO  10  1=1, NN 
pn  in  J= 1  *  NN 
10  B  (  I  » J  )  =  A  (  1  *  J  ) 

INITIALIZE  DETERMINANT  VALUE  AND  ROW  INTERCHANGE  COUNT 

D=1 .0 
KCT  =  0 

PFRFORM  ELIMINATION  ON  N  COLUMNS 
DO  90  1=1, NN 

SEARCH  I-TH  SUB-COLUMN  FOR  I-TH  PIVOT  ELEMENT 
TR^=0. 

DO  BO  T I = I »NN 

TP  ( ARS f R (  T T  ,  T1  1 -TPP )  30,20,20 

?n  MPR=IJ 

tpf=abs(B( it , n  ) 

30  CONTINUE 

IF  PIVOT  ELEMENT  IS  ZFRO,  THEN  DET(A,N)=0.0 

IF  (FUNPR.I))  33,32,35 
32  0=0. 

PO  TO  qs; 

DIVIDP  PIVOT  ROW  BY  PIVOT  FLFMFNT 

3^  no  40  J  =  T , NN 
40  °TVtJT=B(NPP,J)/B(NPR,T ) 

UPDATE  THE  PRODUCT  OF  PIVC"  ELEMENTS  AND  SUM  OF  ROW  INTERCHANGES 


DE  T .0001 
DET . 0002 
DE  T .0003 
DET .0004 
DET. 0005 
DET. 0006 
DET. 0007 
DET. 0008 
DET. 0009 
DET. 0010 
DET. 0011 
DFT. 0012 
DET. 0013 
DET. 0014 
DET. 0015 
DET. 0016 
DET. 0017 
DFT. 0018 
DET. 0019 
OFT. 0020 
DET. 0021 
DET. 0022 
DET. 0023 
DET. 0024 
DET. 0025 
DET. 0026 
DET. 0027 
DET. 0028 
DFT. 0029 
DET. 0030 
DFT. 0031 
DET. 0032 
DET. 0033 
DET. 0034 
DET. 0035 
OFT. 0036 
DET. 0037 
DET. 0038 
DET. 0039 
DET. 0040 
DET. 0041 
DET. 0042 
DFT. 0043 
DET. 0044 
DFT. 0045 
DET. 0046 
DET. 0047 
DET. 0048 
DET. 0049 
DET. 0050 
DET. 0051 
DET. 0052 
DET, 0053 
DFT. 0054 
DET. 0055 
DET. 0056 
DET. 0057 
DET. 0058 
DFT. 0059 
DET. 0060 
DET. 0061 
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D  =  D*B ( NPR  *  I  ) 

KCT=KCT+(NPR-I ) 
r 

r  ELIMINATE  REMAINING  ELEMFNTS  IN  I-TH  SUB-COLUMN 

r 

KTO=NN 
v rooM=NM 
DO  90  K=I *NN 
IF  (KFROM-NPR)  70.80.70 

70  RLE=-B< KFROM, I ) 

DO  75  J= I  *  NN 

75  B(KTO»J)=B(KFROM»J)+RLE#PIV(J) 

KTO=KTO-l 

80  KFROM=KFROM-l 

90  CONTINUE 
r 

r  IF  TOTAL  NO.  OF  ROW  INTFRCHANGES  WAS  ODD.  THEN 

r  NFGATF  THF  PRODUCT  OF  THE  PIVOT  FLEMFNTS 

r 

os  IF  (  KCT  .NE.  2* ( KC  T/2 )  )  D  =  -D 

DET  =  D 
RETURN 
C 

C  GIVE  ERROR  MESSAGE  FOR  INCORRECT  VALUE  OF  N 

r  AND  RETURN  TO  SYSTEM  VIA  EXEM 

r 

IPO  WRITE  (6.1000)  MN 
CALL  FXEM 
rctuRN 

r 

1000  FORMAT  (3H0N=. I 12.80H  IS  INCORRECT  FOR  FUNCTION  DET) 
Fr\|n 


DET. 0062 
DET. 0063 
DET. 0064 
DET. 0065 
DET. 0066 
DET. 0067 
DET. 0068 
DET. 0069 
DET. 0070 
DET. 0071 
DET. 0072 
DET. 0073 
DET. 0074 
DET. 0075 
DET. 0076 
DET. 0077 
DET. 0078 
DET. 0079 
DET. 0080 
DET. 0081 
DET. 0082 
DET. 0083 
DET. 0084 
DET. 0085 
DET. 0086 
DET. 0087 
DET. 0088 
DET. 0089 
DET. 0090 
DET. 0091 
DET. 0092 
DET. 0093 
DET. 0094 
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SIBFTC  MSUB  DECK 

CMSUBJ  VERSION  1  OF  RS  MSUB 

SUBROUTINE  SIMPLX  ( I NF I  X , A , B , TOL * PRM , KOUT * ERR . JH , X *P » Y , KB , E ) 

C 

DIMENSION  INFIX! 8) ,A (1 ) ,6(1 ! .TOL (4) »KOUT (7 )  »ERR ( 8 )  , JH ( 1  )  »X ( 1  )  * 

1  P( 1) »Y( 1 ) ,KB( 1 ) *F( 1 ) ,ZZ(3) »  I  OF  I  X ( 16 )  *  TFRRC8) 

C 

FOU I  VALENCE  ( I  NFL AG , I  OF  I X i 1 )  )*  (N  ,  I  OF  I X ( 2 )  !  , 

1  ( ME  » IOF I X ( 3 )  ),  (  M » I  OF  I  X ( 4 )  )»  < MF , I  OF  I X ( 5 >) , 

2  t  MC  *  IOF IX! 6!  )»  (  NCU1,  IOF IX! 7)  1  ,  l  NvER  ,  I  OF  I x t 8 5  ), 

3  {  K  *  I  OF  1X19!  !  ,  (IlER,  IOFIxilu)  ),  (INvC  ,  lOFIxill)  )  , 

4  ( NUMVR ,  IOF IX 112)  ),  l  NuMPv,  IOFIa(13)  )  , 

5  (INFS,  1  OF  I X  (  1 4 )  )  ,  !  JT,  I  OF  I X  (  15  J  )  ,  (  LA  ,  IOFIXU6)  ), 

6  ( ZZt 1 ) »TP IV )  ,  (ZZ< 2) ,TZERO) , l ZZ( 3 ) ,TCOS I  ) 

MOVE  INPUTS  ...  ZERO  OUTPUTS 

DO  1340  1=  1,  8 

TERR ! I )  =  0.0 
I  OF  I  X ( I +8  )  =  0 

1340  I  OF  I  X  t  I  )  =  INFIXCI) 

LA  =  0 

DO  1308  1=1,3 

1308  ZZ(  I  )  =  TOL! I ) 

TCOST  =  -  ABS  ( T  COST ) 

PMI X  =  PRM 
M2  =  M**2 
INFS  =  1 

CHECK  FOR  ILLEGAL  INPUT 
IF  (N)  1304,  1304,  1371 

1371  IF  (M  -  MF  )  1304,  1304,  1372 

1372  IF  (MF  -  MC)  1304,  1304,  1373 

1373  IF  (  MC  )  1304  ,  1304,  1374 

1374  IF  (ME  -  M  )  1304,  1375,  1375 

1375  IF!  MOD  (INFLAG,  4  )  -  1  )  1400,  1320,  100 


NEW  1 

*****  SUBROUT 


STARTS  PHASE  C  IE 

INE  NEW  (M,N ,  JH,  KB,  A,  B,  MF ,  ME  ) 


INITIATE 

1400  DO  1401  I  =  1,  M 
1 4nl  JH ( 1 )  =  0 

C  INSTALL  SINGLETONS 

KT  =  0 

DO  1402  J  =  1,  N 
KB ( J )  =  0 
MM  =  KT  +  MF 
LL  =  KT  +  M 

C  TALLY  ENTRIES  IN  CONSTRAINTS 

KQ  =  0 

DO  1403  L  =  MM  ,  LL 

IF  t A ! L  »  J  1404,  1403,  1404 

1404  KO  =  KQ+1 
LO  =  L 

1403  CONTINUE 

C  CHECK  WHETHER  J  IS  CANDIDATE 

IF  (KO  -  1)  1402,  1405,  1402 

1405  TA  =  LO-  KT 

IF  !  JH(IA)  )  1402,  1406,  1402 

1406  IF  ( A  I  LQ ) *B ( I A ) )  1402,  1407,  1407 

C  J  IS  CANDIDATE.  INSTALL 


MSUB. 001 
MSUB. 002 
MSUB. 003 
MSUB. 004 
MSUB. 005 
MSUB. 006 
MSUB. 007 
MSUB. 008 
MSUB. 009 
MSUB. 010 
MsuB.Oll 
MouB.u 1 2 
MSUB. 013 
MSUB. 014 
MSUB. 015 
MSUB. 016 
MSUB. 017 
MSUB. 018 
MSUB. 019 
MSUB. 020 
MSUB. 021 
MSUB. 022 
MSUB. 023 
MSUB. 024 
MSUB. 025 
MSUB. 026 
MSUB. 027 
MSUB. 028 
MSUB. 029 
MSUB. 030 
MSUB. 031 
MSUB. 032 
MSUB. 033 
MSUB. 034 
MSUB. 035 
MSUB. 036 
MSUB. 037 
MSUB. 038 
MSUB. 039 
MSUB. 040 
MSUB. 041 
MSUB. 042 
MSUB. 043 
MSUB. 044 
MSUB. 045 
MSUB. 046 
MSUB. 047 
MSUB. 048 
MSUB. 049 
MSUB. 050 
MSUB. 051 
MSUB. 052 
MSUB. 053 
MSUB. 054 
MSUB. 055 
MSUB. 056 
MSUB. 057 
MSUB. 058 
MSUB. 059 
MSUB. 060 
MSUB.Otsl 


46 


AFML-TR-67-121 
PART  III 


1407  JH ( IA)  =  J 
KB(J)  =  I A 
3402  KT  =  KT  +  ME 
r 

r**^ND  OF  NPW 

r 

r 

1320  CONTINUF 

r 

C  VER  1  FORMS  INVERSE  FROM  KB 

C*****SUBROUTINE  VER  (  A,  B,  JH,  X,  E,  KB*  Y*  IOFIX.  TPIV*  M2 
C 

C  INITIATE 

1100  ASSIGN  110?  TO  KP I V 
ASSIGN  1114  TO  KJMY 

IF  (LA)  1121,  1121,  1122 

1121  I N VC  =  0 

1122  NUMVR  =  NUMVR  +1 

DO  1101  I  =  1,  M2 

1101  E( I  1=0. 

MM=  1 

DO  1113  1=1,  M 

F ( MM )  =1.0 
X  (  I  )  =  B  (  I  ) 

1113  MM  =  MM  +  M  +  1 

DO  1110  I  =  MF,  M 

if  (  jh  ( in  mi,  mo,  ini 

1111  JH ( I )  =  12345 
1110  CONTINUF 
INFS  =  1 

FORM  INVERSE 
oo  no?  JT=  1  ,  N 
IF  (  KB(JT))  600  ,  1102  ,  600 
C  600  CALL  JMY  ( JT ,  A,  E,  M,  Y  ) 

C  CHOOSE  PIVOT 

1114  TY  =  0. 

DO  1104  I  =  MF,  M 

IF  ( JH  < I  )  -  12345  )  1104,  1105,  1104 

1105  IF  (  ABS  (  Y ( I )  )  -  TY  )  1104,  1104,  1106 

1106  I R  =  I 

TY  =  ABS  (  Y ( I  )  ) 

1104  CONTINUE 

C  TEST  PIVOT 

IF  (TY  -  TPIV  )  1107,  1108,  1108 

C  BAD  PIVOT,  ROW  IR,  COLUMN  JT 

1107  KB ( JT ) =  0 
GO  TO  1102 

C  PIVOT 

1108  JH< IR)  =  JT 
KB(JT)  =  IR 
GO  TO  900 

C  900  CALL  PIV  (IR,  Y,  M,  E,  Z,  X) 

1102  CONTINUE 

r  RESET  ARTIFICIALS 

DO  1109  I  =  1,  M 

IF  (  JH ( I )  -  12345  )  1109,  Ill?,  1109 

n 12  jh( n  =  o 

1109  -CONTINUE 
OF  VFR 

C 

C 


MSUB.062 
MSUB.063 
MSUB.064 
MSUB.065 
MSUB.066 
MSUB.067 
MSUB.068 
MSUB.069 
MSUB.070 
MSUB.071 
MSUB.072 
MSUB.073 
MSUB.074 
MSUB.075 
MSUB.076 
MSUB.077 
MSUB.078 
MSUB.079 
MSUB.080 
MSUB.081 
MSUB.082 
MSUB.083 
MSUB.084 
MSUB.085 
MSU3.086 
MSUB.087 
MSUB.088 
MSUB.089 
MSUB.090 
MSUB.091 
MSUB.092 
MSUB.093 
MSUB.094 
MSUB.095 
MSUB.096 
MSUB.097 
MSUB.098 
MSUB.099 
MSUB.100 
MSUB.101 
MSUB. 102 
MSUB.103 
MSUB. 104 
MSUB. 105 
MSUB. 106 
MSUB. 107 
MSUB. 108 
MSUB. 109 
MSUB. 110 
MSUB. Ill 
MSUB. 112 
MSUB. 113 
MSUB. 114 
MSUB. 115 
MSUB. 116 
MSUB. 117 
MSUB. 118 
MSUB. 119 
MSUB. 120 
MSUB. 121 
MSUB. 122 
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100 

ASSIGN 

705 

TO  NOEL 

MSUB. 123 

ASSIGN 

1000 

TO  KJMY 

MSUB. 124 

ASSIGN 

221 

TO  KPIV 

MSUB. 125 

r 

MSUB. 126 

C 

PERFORM  ONE  ITERATION 

MSUB, 127 

r 

MSUB. 128 

c  xck 

:  i  x 

CHECKER 

MSUB. 129 

C*****SU8ROUT INE  XCK  (  M,  MF , 

JH,  X,  T2ERO,  JIN  ! 

MSUb. 130 

C 

MSUB. 131 

C 

RESET 

X  AND  CHECK 

FOR  INFEASIBILI T IES 

MSUB. 132 

1200 

JIN  = 

0 

MSUb. 133 

NEG  =  0 

MSUB. 134 

DO  1201 

I  = 

MF,  M 

MSUB. 135 

IF  <  ABS 

( 

XII)  )  -  T ZERO )  1202,  1203,  1203 

MSUB.  136 

1202 

X(I)  =  0 

.0 

MSUB. 137 

GO  TO  1201 

MSUB. 138 

1  203 

IF  (  XII)  ) 

1208,  1201, 

1205 

MSUB. 139 

1205 

IF  (  JH ( I )  ) 

1201,  1206 

,  1201 

MSUB. 140 

1208 

NEG  =  1 

MSUB. 141 

I  206 

JIN  = 

1 

MSUB. 142 

1201  CONTINUE 
C**END  OF  XCK 
C 

r 

C  CHECK  CHANGE  OF  PHASE..  GO  BACK  TO  INVERT  IF  GONE  INFEAS. 

IF  (INFS  -  JIN  )  1320,  5C0,  200 

C  BECOME  FEASIBLE 

200  INFS  =  0 
2 Cl  PM  I  X  =  0.0 

r 

C  GET  1  GET  PRICES 

C*****SUBROUTINE  GET  (  M,  MC»  MF,  JH,  X,  P,  E,  INFS,  PMIX  ) 
r 

500  MM  =  MC 

C  PRIMAL  PRICES 

502  DO  503  J  =  1,  M 
P { J  >  =  F ( MM  ) 

503  MM  =MM  +  M 

IF  {  INFS  )  501,  599,  501 

r  COMPOS  I TF  PRICES 

501  no  504  J  =  1  ,  M 

504  PfJJ  =  P{J)*  PMIX 
DO  505  I  =  MF,  M 

MM  =  j 

IF  t  XII)  )  506,  507,  507 

506  DO  508  J  =  1,  M 
P(J)  =  P ( J !  +  El  MM  ) 

508  mm  =  MM  +  M 
GO  TO  505 

507  ;f  1 JH ( I) )  505,  509,  505 

509  DO  51A  J  =  1 ,  M 

P  (  J  )  =  P  C  J  )  -  E  (  MM  ) 

519  MM  =  MM  +M 

505  CONTINUE 
C 

599  CONTINUE 
C**END  OF  GET 
r 


MIN  D-J.  SELECTS  COLUMN  TO  EN1ER  BASIS 
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MSUB. 160 
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MSUB, 162 
MSUB. 163 
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*****  SUBROUT  I  ME  MIN  (  JT  »  N»  M»  A,  P,  KB,  ME,  I  COi>  I  ) 

700  JT  =  0 

BB  =  T  COST 

701  DO  702  JM  =  1,  N 

SKIP  COLUMNS  IN  BASIS 
703  IF  (  KB(JM)  )  702,  300,  702 

300  CALL  DEL  (  JM,  DT,  M,  A,  P) 

705  IF  (  DT  -  BB  )  708,  702,  702 

708  BB  =  DT 
JT  =  JM 

702  CONTINUE 


*END  OF  MIN 


203 


IF  (JT)  203,  203,  600 

ALL  COSTS  NON-NEGATIVE...  K  = 
K  =  3  +  INFS 
60  TO  257 

NORMAL  CYCLE 


3  OR  4 


r  JMY  l  J  MULTIPLY.  BASIS  INVERSE  *  COLUMN  JT 

C*****SUBROUTINE  JMY  ( JT ,  A,  E,  M,  Y,  ME  ) 
r 

600  DO  610  I =  ] ,M 
610  Y(I)  =0. 

I_d  =  JT *ME  -  ME 
LL  =  0 

DO  605  I =  1 ,M 
LP  =  LP  +  1 

IF  ( A  C  LP  )  )  601,  602,  601 

601  DO  606  J  =  1 , M 

LL  =  LL  +  1 

606  Y ( J )  =  Y ( J )  +  A ( LP )  *  E(LL) 

GO  TO  605 

602  LL  =  LL  +  M 
605  CONTINUE 


699  GO  TO  K JMY  ,  (  1000  ,  1114  ,  1392  ) 

'"**END  OF  JMY 

r 

r 

C  ROW  1  ROW  SELECTION — COMPOSITE 

0***** SUBROUT  I NE  ROW  (  IR,  M,  MF ,  JH,  X,  Y,  I P I V  ) 
r 

C  AMONG  EQS •  WITH  X=0,  FIND  MAX  ABS(Y)  AMONG  ARTIFICIALS,  OR,  IF  NONE, 
r  GET  MAX  POSITIVE  Y (  I  )  AMONG  REALS. 

1000  IR  =  0 

A  A  =  r> .  f) 

I A  =  0 

DO  1050  I  =  mf,  m 

IF  (  X(II  )  1050,  1041,  1050 

1041  Y I  =  ABS  (  Y (  I  )  ) 

IF  (  Y I  —  TPIV  )  1050,  1050,  1042 

1042  IF  (  JH ( I )  )  1043,  1044,  1043 

1043  IF  (14)  1050,  1048,  1050 

1048  IF  (  Y(I)  )  1050,  1050,  1045 

1044  IF  (IA)  1045,  1046,  1045 

1045  IF  (  Y I  -  AA  )  1050,  1050,  1047 

1046  I A  =  I 
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MSUB. 194 
MSUB. 195 
MSUB. 196 
MSUB. 197 
MSUB. 198 
MSUB. 199 
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MSUB. 203 
MSUB. 204 
MSUB. 205 
MSUB. 206 
MSUB. 207 
MSUB. 208 
MSUB. 209 
MSUB. 210 
MSUB. 211 
MSUB. 212 
MSUB. 213 
MSUB. 214 
MSUB. 215 
MSUB. 216 
MSUB. 217 
MSUB. 218 
MSUB. 219 
MSUB. 220 
MSUB. 221 
MSUB. 222 
MSUB. 223 
MSUB. 224 
MSUB. 225 
MSUB. 226 
MSUB. 227 
MSUB. 228 
MSUB. 229 
MSUB. 230 
MSUB. 231 
MSUB. 232 
MSUB. 233 
MSUB. 234 
MSUB. 235 
MSUB. 236 
MSUB. 237 
MSUB. 238 
MSUB. 239 
MSUB. 240 
MSUB. 241 
MSUB. 242 
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1047  4 A  =  Y I 
IR  =  I 
1050  CONTINUE 

IF  (  IR)  1099,1001 .1099 

1001  AA  =  1.0E+20 

C  FIND  MIN.  PIVOT  AMONG  POSITIVE  EQUATIONS 

DO  1010  IT  =  MF  ,  M 

IF  C  Y  <  I  T  )  -  TP  I  V  )  1010*  1010*  1002 

1002  IF  <  X ( I T  )  )  1010,  1010*  1003 

1003  XY  =  X< IT )  /  YCIT) 

IF  (  XY  -  AA  )  1004*  1005,  1010 

1005  IF  <  JH ( I T ) )  101C,  1004,  1010 

1004  aa  =  XY 
IR  =  IT 

3010  CONTINUE 

IF  (NEG)  1016,  1099,  1016 

C  FIND  PIVOT  AMONG  NEGATIVE  EQUATIONS,  IN  WHICH  X/Y  IS  LESS  THAN  THE 
C  MINIMUM  X/Y  IN  THE  POSITIVE  EQUATIONS,  THAT  HAS  THE  LARGEST  AdSF(Y) 
1016  BB  =  -  TPIV 

DO  1030  I  =  MF  ,  M 
IF  (XIIII  1012,  1030,  1030 
1012  IF  <  Y ( I  )  -  BB  )  1022,  1030,  1030 

1022  IF  (  Y  (  I  )  *  AA  -  X(I)  )  1024,  1024,  1030 

3  024  BB  =  Y (  I  ) 

IR  =  I 

1030  CONTINUE 
1099  CONTINUE 
C**END  OF  ROW 

TEST  PIVOT 

206  I  Ft  IR  )  207,  207,  210 

C  NO  PIVOT 

207  K  =  5 

25^  if  ( PM I X )  201,  400,  201 

C  ITERATION  LIMIT  FOR  CUT  OFF 

210  I F  ( I TFR  -NCUT  )  900,  160,  160 

r  PIVOT  FOUND 

r 

C  PIV  1  PIVOT.  PIVOTS  ON  GIVEN  ROW 

C*****SUBROUTINE  PIV  (  IR,  Y,  M,  E,  X,  NUMPV,  TECOL  ) 

C  LEAVE  TRANSFORMED  COLUMN  IN  Y(I) 

C 

900  NUMPV  =  NUMPV  +  1 

Y I  =  -Y(IR) 

Y ! IR)  =  -1 . 

LL  =  0 

r  TRANSFORM  INVFRSF 

903  DO  904  L  =  IR,  M2,  M 

IF  C  F(L)  )  905,  914,  905 

914  LL  =  LL  +  M 
GO  TO  904 

905  XY  =  ECL)  /  Y I 
F  t  L  !  =0. 

DO  906  I  =  1 ,  M 
LL=  LL  +1 

906  E  <  LL )  =  ECLL)  +XY*  Y!IJ 

904  CONTINUE 

r  TRANSFORM  X 

XY  =  X( IR)  /  VI 
X ( IR)  =  0. 

DO  908  I  =  1,  M 
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908  X( I)  =  X( I)  +XY*  Y ( I ) 

C  RESTORE  Y< IR) 

Y( IR)  =  -Y I 
C 

999  GO  TO  KPIV  ,  (  221.  1102  ) 

*END  OF  PIV 

221  IA  =  JH(IR) 

IF  (  IA  )  213,  213,  214 

214  KB|  !A  )  =  0 

213  KB(JT)  =  IR 

JH( IR)  =  JT 

LA  =  0 

ITER  =  ITER  +1 

I NVC  =  I NVC  +1 

INVERSION  FREQUENCY 
IF  ( 1 NVC  -  NVER  )  1200,  1320,1200 

CUT  OFF  ...  TOO  MANY  ITERATIONS 

160  K  =  6 


ERR  1 

*****SUBROUT 


ERROR  CHECK.  COMPARES  AX  WITH  B,  PA  WITH  ZERO 
INE  ERR  (  M,  A,  B,  TERR,  JH ,  X,  P»  Y,  ME,  LA  ) 


STORE  AX-B  AT  Y 

400  ASSIGN  410  TO  NDEL 
DO  401  I  =  1,  M 

401  Y( I)  =-B( I  I 

DO  402  I  =  1,  M 
JA  =  JH ( I ) 

IF  (JA)  403,  402,  403 
403  IA  =ME*  ( J A— 1 ) 

DO  405  IT  =  1,  M 
I A  =  I A  +  1 

I F  <  A ( I A )  )  415,  405,  415 

415  Y(  IT)  =Y ( I T)  +X( I )  *  A ( I  A ) 

405  CONTINUE 

402  CONTINUE 

FIND  SUM  AND  MAXIMUM  OF  ERRORS 

DO  481  I  =  1,  M 
Y I  =  Y  (  I  ) 

IF  (  JH  < I )  )  472,  471,  472 

471  Y I  =  Y  I  +  X  (  I  1 

472  TERR ( LA  +  1 )  =  TERR ( LA+1 )  +  ABS  (YI) 

IF  (  ABS  ( TERR ( LA+2 ) ) -  ABS  (  YI  )  )  482,  481,  481 

482  TFRR (LA+2 )  =  YI 
481  CONTINUE 

STORE  P  TIMES  BASIS  AT  DT 

DO  411  I  =  1,  M 
JM  =  JH ( I ) 

IF  (  JM  )  300  ,  411  ,  300 

300  CALL  DEL  (  JM,  DT,  M,  A,  P) 

410  TERR ( LA+3 )  =  TERR ( LA  +3)  +  ABS  ( DT ) 

IF  (ABS  ( TERR ( LA+4 ) )  -  ABS  ( DT )  )  413,  411,  411 

413  TERR (LA+4 )  =  DT 

411  CONTINUE 
*END  OF  ERR 


IF  (LA)  193,  191,  193 

191  LA  =  4 
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IF  (INFLAG  -  4  1  1320.  193.  193 

193  IF  (K— 5)  1392.  194.  1392 

194  ASSIGN  1392  TO  KJMY 
GO  TO  600 

600  CALL  JMY  (.....) 

GO  TO  1392 
1304  K  =  7 

SFT  FXIT  VALUES 
1  3Q7  DO  1309  1=  1,  8 

1  309  FRR  (  I  )  =  TERR ( I ) 

00  1329  1=1.7 

1329  KOUT ( I )  =  I  OF  I X ( 1+8) 

RETURN 


DEL 

****  *SUBROUT 


DELTA-JAY.  PRICES  OUT  ONE  MATRIX  COLUMN 
INE  DEL  (  JM.  DT,  M.  A.  P.  ME  ) 


300  DT  =  0. 

LL  =  (JM  -  1)  *  ME 

301  DO  303  MM  =  1.  M 
LL  =  LL  +  1 

IF  (  A(  LL  ) ) 304 .  303.  304 
304  DT  =  DT  +  Pt  MM  )  *  A  (  LL  > 
303  CONTINUE 

399  GO  TO  NDEL  .  (  410  .  705  ) 
**END  OF  DEL 

END 
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Figure  1.  Resulting  MWD  Using  Ten  Sets  of  Molecular  Weights  a 
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Sets  of  Molecular  Weights  and  g  =  4. 0,  (Compare  with  Figures  1,  3,  and  4) 
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Figure  8.  The  Effect  of  Using  More  Experimental  Data  on  the  MWD 
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